{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" }, "colab": { "name": "804_Heat Equation-CN-Different_Initial_Boundary.ipynb", "provenance": [], "include_colab_link": true } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "0uTl1O_FDeba" }, "source": [ "# The Implicit Crank-Nicolson Difference Equation for the Heat Equation\n", "## The Heat Equation\n", "The Heat Equation is the first order in time ($t$) and second order in space ($x$) Partial Differential Equation: \n", "\\begin{equation} \\frac{\\partial u}{\\partial t} = \\frac{\\partial^2 u}{\\partial x^2},\\end{equation}\n", "The equation describes heat transfer on a domain\n", "\\begin{equation}\\Omega = \\{ t \\geq, 0\\leq x \\leq 1\\}. \\end{equation}\n", "with an initial condition at time $t=0$ for all $x$ and boundary condition on the left ($x=0$) and right side ($x=1$).\n", "\n", "## Crank-Nicolson Difference method\n", "This note book will illustrate the Crank-Nicolson Difference method for the Heat Equation with the __initial conditions__ \n", "\\begin{equation} u(x,0)=x^2, \\ \\ 0 \\leq x \\leq 1, \\end{equation}\n", "and __boundary condition__\n", "\\begin{equation} u(0,t)=t, u(1,t)=2-exp(-t). \\end{equation}\n" ] }, { "cell_type": "code", "metadata": { "id": "JI8r-JYhDebd" }, "source": [ "# LIBRARY\n", "# vector manipulation\n", "import numpy as np\n", "# math functions\n", "import math \n", "\n", "# THIS IS FOR PLOTTING\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")\n" ], "execution_count": 1, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "kp77vORCDebf" }, "source": [ "## Discete Grid\n", "The region $\\Omega$ is discretised into a uniform mesh $\\Omega_h$. In the space $x$ direction into $N$ steps giving a stepsize of\n", "\\begin{equation} h=\\frac{1-0}{N},\\end{equation}\n", "resulting in \n", "\\begin{equation}x[i]=0+ih,\\end{equation}\n", "and into $N_t$ steps in the time $t$ direction giving a stepsize of \n", "\\begin{equation}k=\\frac{1-0}{N}\\end{equation}\n", "resulting in \n", "\\begin{equation}t[i]=0+ik.\\end{equation}\n", "The Figure below shows the discrete grid points for $N=10$ and $N_t=15$, the red dots are the unknown values, the green dots are the known boundary conditions and the blue dots are the known initial conditions of the Heat Equation." ] }, { "cell_type": "code", "metadata": { "id": "aUpklLOODebf", "outputId": "c4aee95e-81c6-4644-91da-3d8bbc973973", "colab": { "base_uri": "https://localhost:8080/", "height": 441 } }, "source": [ "N=10\n", "Nt=100\n", "h=1/N\n", "k=1/Nt\n", "r=k/(h*h)\n", "time_steps=50\n", "time=np.arange(0,(time_steps+.5)*k,k)\n", "x=np.arange(0,1.0001,h)\n", "X, Y = np.meshgrid(x, time)\n", "fig = plt.figure(figsize=(8,6))\n", "plt.plot(X,Y,'ro');\n", "plt.plot(x,0*x,'bo',label='Initial Condition');\n", "plt.plot(np.ones(time_steps+1),time,'go',label='Boundary Condition');\n", "plt.plot(x,0*x,'bo');\n", "plt.plot(0*time,time,'go');\n", "plt.xlim((-0.02,1.02))\n", "plt.xlabel('x')\n", "plt.ylabel('time (ms)')\n", "plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n", "plt.title(r'Discrete Grid $\\Omega_h,$ h= %s, k=%s'%(h,k),fontsize=24,y=1.08)\n", "plt.show();" ], "execution_count": 2, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn0AAAGoCAYAAAAzaFPMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfZzUdb3//+drd1lgBfECUq4XhBV2RURI5fTz6mhFncTKVIo0NUWXYyWd7KCWmmZZ1rFD3xTFNC+20jwnQ+t0ZWlZWkEIynIRGhdykQiKIpfLvn5/fD4Ts8vM7MzuvOa9r88877fb+wY7Ozv7YBb9vJn5vGZEVUFEREREyVYROoCIiIiI7HHTR0RERFQGuOkjIiIiKgPc9BERERGVAW76iIiIiMoAN31EREREZYCbPiIiIqIywE1fCYjIUyKiInJR6BYqjq7+TOOvVRGpLWqYEyKyOv7znxa6pRzx/icqT9z05SAi30s7OKfWXhHZIiKrROQxEblWREaEbg1NRC4SkRtF5Lhu0NJLRC4Tkf8Rkb+LyHYR2S0im0TkNyJys4gcG7qzq0RknIj8t4i8KCJviMhOEVkjIj8SkfNEREI3JomIHBnf3y+JyC4R+YeIPC4iZ3Ty9vqKyNT47+P/ichraf+fGVPs/u6s2PdtV2+bPxtKLFXlyrIAfA+AAtgDYFO8/gFgZ3x5arUCeARA/yy38wCA5QA+FPrPZHhfPRXfFxcF7jgLwMZ2P5+dALYC2Nfu8l9k+5nl8X269DNNa6jtxNdWAZjT7s+zE8C2dn++hQBGhv67keXPsDpuPC10S569xwJ4Le2+3ZZ2/7cCmN2J2/xgu59X+hpTLve/xX3b1dsO+bPh4rJcfKQvP39U1SPjdYSq9gZwKID3AXgY0f8IzgXwvIgMaf/Fqnqhqo5R1R+XNru8iMgnATwG4EgAKwBcAmCQqvZW1cMAVAOYCOCLADYAeA+AA35e+Qj1M40fvfsRgE8hOmjdDuDo+M/YD8A74s+9BuB4AM+JSF0pG5NGRHoDmA/gcACLABwT39eHAvgmAAHwFRF5Tydu/lUAPwPwJQAzilPsh+V9W4TbLuufDSVU6F1nd17Y/0jfUx1cbwr2P/r3XOjuQPfVUwj4SB+iDc7uuOExAL06uH41gJsBNATq7dQjfQCuwf5Hn9+X43q12P9ozosAqkP/HWnXl2o7LXRLHq1Xxa1vARic4fM/jj+/sMDbrczwMyurR/qs7tuu3nbInw0Xl+XiI31FoKo/B/C5+MMTReSs9M/nOulfRKpF5DMi8sf4vKy98Tkni0XkOyIyOdP3FJGxIjJXRFaKyI74a18QkTkiMrHddf950raIDBaRO0Tk5fg8t+cz3PYxInJvfD7crvi2/yAiV4hIj3bXvUhEFMCp8UX3tTsHcnVXbr8AX0a0kVsD4OOquivXlVV1j6p+UVWXtmvL677K9TONP18hIp+Kf447RWRzfB5Rxp9nPkTkUESbPgC4SVX/L8efbzWA6fGHDQAydnYHInKYiPxX/Pdht4isF5F5IjIwdFssdT9+X1XXZ/j8bfGvx4vI0fneqKru63KZARE5VESejf9+LxaRIwy/ncl929Xb7q4/G6IuC73r7M4LeT7SF1+3GtH5fgrgB+0+9xQyPAqG6Nys1OdS55i8DqAl7bIfZvhen2p3ne3x12mmXuz/V/0MAJvj378df93z7a57JdqeK/ZWu+/1WwA1adc/H9G5jnuw/5yZTWnrL125/Tx/TkPTvn5WF3/med1X2X6maT/Xx9Ka9qb9fPYC+HDa52oLaLsi/po3AfTN82t+G3/NX/L9PiX6byt1P3887fdvA9iVdt/8HcChgTv7xv9dKoAPZ7lOBYA34uvM7ML3qk37swd5pA/RqRFL4s89a3n/W963xb7tUv5suLgsFx/pKxJV3QPgyfjDk/P8so8heoRsB4ALEG12DgXQE8BwRBukxelfICLnIjqJvxLAowDqVbVP/HWHIzqILszy/b6JaMjhXap6kKr2AfCRtNv+IIBvIzr4fh7AAFXtC6AG0VPYfwNwGqLzyFJ/7odV9UgAf4wv+ozuP//xSFV9Z1duP0+npv3+pwV+bTY576sO/CeAsxEddK4G0C/++YwE8GsA93ay6fT419+r6lt5fs3P4l+PF5F+nfy+lr6NaEP8L6p6EIA+iO67NxAdaK/J/qUlMRbRuV8AsDTTFVS1FdE5pABQX4ooCyIyHMDvAYxD9P+yM1X1dcNvaXnfls3PjagQ3PQV1wvxr4PzfJrypPjXB1T1IY2fklTVfaq6VlW/o6pfTV05vs3UhugHqnquqi5LfV5Vt6pqk6r+R5bv1wLg3ar6x7SvWRXfdiWAb8UXn6uqt6nqa/F19qjqLxANruwAcEmhT70Z3/7Y+NddiDaOxZD1vspFRA5CtOkDgJtV9RuquiP++r8jmgrM9FRTPkbGv67Iea22UvdHBaJHRLub3Yg2F88CgKq2qOp8RE/XAxk22hK9NFD7l1LKd91YYF/638MNOa6X+lx3eUq6IPHLkDwDYBSAnwD4N1V9O8P1innfW963ZfFzIyoUN33Flf6v4sPyuP6b8a/5/g/nDACDET09enUBXSkPqOo/snzuNESPLr4Yb8AOoKovAXgO0dOXpxX4vS1vP3Vfv6GqmukKInKdRK/T1379d5bbzHVf5fIeRE8t7UaGRyxVdTeAb3TidgEg9UjdAQfjHLan/f7g1G9EpLeI7BOR6zvZUix3q+qWDJc/Fv86It5Ip9uO6FSKzqztKEz6996Z43o74l/7FHj7wYnI8Yge4RsCoAnAR+K/p5kU8763vG8T/3Mj6oyq0AFl7v8QPxUoIvMRnUP4dJaDILD/kcHFmvnE5I48m+Nz/xL/OlpENuW4XmrjUeijRta335G+ADKdkJ7tKc9c91Uux8e/Pq+q27Jc5+lO3vYb8a+HFPA16dd9Le33DYj+0fdiJ1uK5S9ZLk//+30I0ja6qvoNdH7jTG2djGgI7WAAdwL492z/cAJ43xN5x0f6iuvQtN9v7ejKqvo0gOsRPZV4FoD/AfCaiCwTkW+IyOh2X5LatKztZN/mHJ9LPdrYM/4+2Vav+Ho1BX5vy9tP3deHiGR+FwpVna2qklqIHtHIJdd9lcuA+NdcTyl19undV+JfGwr4mtR1FW2bxsW/ht70ZTw3UdtOX3d2orsY0h9V7Z3jeqm/r4U+khjaTYg2fL9R1Zm5NnwGLO/bpP/ciDqFm77iSh1IX1HVvfl8gareDKAO0Qnrv0D0lO8YAP8BoFlELixiX66XIUj9XfhJ+uYox7qxwO9tefup8xp7AWi/Ue6s7viSDb+Of32XiOT7aN+/xb8+q6rpB7ZxiM6B7PA8xTKXvlEelON6qc9tNGyx8HD867+KSGOJv7flfZv0nxtRp3DTVyQiUo3onDsgOj8mb6r6d1W9VVWnIDo/7XQAv0P09PsdIvKO+Kqpc8yGFyG5vdRtDzO4bevbT3+69N+yXqs0Uo8Q5nOgKdT/InppnGoAHR6g49cETE1P/6Ddp8ch2iw3iMh8EdkWn+N4XSfbSkZEPpfl/Mx81uc6/g5tLEf0KCmQ5RFWEakAkHqdt+ZO/rFCmQvgs/HvvyMil+S6cpHve8v7Nuk/N6JO4aaveC5D9BZYQMdPHWYVT+4+BeADiF7T7SAAk+JPPxf/eqyIDO7s98gidQ5bZ2+7Nf4149OrRbj9rFR1HaLzIwHgMyIS8qTsv8a/HiciB2e5zqlZLs9JVTcB+Fr84RdEJOvTvPHww3fjD18CcE+7q4xD9MjoLxEd8D4fX+/LInJKZ/pKqA9ynyKQaxX0dyN+aZwF8YfvznK1E7H/3NAns1yn21LV2wHMRvTf7jwR+XiOqxftvre8b8vh50bUGdz0FYGIvBf7X939WVXN67Xi4kcHs9mD/U8x9ox/fRLR+WCVad+vWJ4EsC6f25bonSHaS00iZ3vasau335EvILrPhgN4SER6dXB9K79EdF/0BPCZ9p+Mf+bZXlInH19B9DqMNQAeFZEDpsTjl8e5B9FL2ewEMC39HDkRGYD9B+J3xec73gVgWnyVSejGVPXGPE8RKMZpCQDw/fjX6VleSij1CNZCVS3k5XS6DVX9GoAbEB0Tvici52W5XrHve8v7NvE/N6JCcdPXSSLST0TeKyI/QPQCuL0RbWryfQFfAHhARO6Lb6dv2m3XArgf0SMxOxE/XRyfJ5jaMHxURB6JX18r9XWHichlIjKn0D9PfNtXInpK5KMi8piIHJd22z1EZJKIfB3ROyW0l3oB1A9LhhcBLsLtd9T/VwAzET3ieDaA50Wkzev9SeQoEbkawHsL/R55drwN4OvxhzeIyGcleuP31M/1x+jCZHK8eXsfotffGwNgkYickPq8iByJaIM9DdE5e+eo6oJ2N3Ns/OuXVPXltMtT56HuaHd9SPS2dKnXWzuts/1O3YXo7f36AnhCROoBQET6xn9fPxxf79r2Xyjx2xTGqzbD5/unFtoOgh2S/rn4qcj2X1vUn4mq3gTgFkT/MGuS6MXUrZndt1257fh6nf7ZEHVb2g3eFqS7Lux/G7Y92P+2Yv9AdFDUtNWK6ITo/llu5ylkfhu2x9rdxuuIps5Sl7UAuCDD7X0WB76VWT5vw3ZaHn/mixG9xlzqtnYA2IK2b5WmGb5uTNrX7UX0iORqAM8U4/YL+Jmdheik7PSfz05E59rtanf5zwAc3Zn7KtvPNP5c0d+GDdFGMf3t7d5Mu42X0653b9rlu9p9zfnxdVJvRD+k3fc4Lb789Azf/7S02+3w71Eef54O7+dC7yPLBWA8ope8STVtw/7/BlsBzM7ydRfl+nO0+/uYa2X62k7/THLd/4geiVdE/53+m9f7tiu33dWfDRdXd138F0p+emD/02GHIzqYvgxgPoDrAIxU1fM1foeJAsxGdC7Vz+Pbq0b0r+yXANwH4HhVfbD9F6nqfwGYEF9nddyniN4z878BzCqwI/2270N0cvO3ED16tw/RSzpsQbTRuQH7T35O/7rliM6d+Tmi/7Eeieip1iHFuP0C+h9H9M4VMxA9qrYG0f/cD0a08Xoa0VOkDar6fjV4WkdVWwCcA+DTiH4mLYj+nD8FcKqq/m8nbrYSbc+P6pv2uYosv2//8jipl64YB+AfqvoK2hof/7okw/dPPWK6A2V40ruqLgZwDKK3QHwZ0X27BdHP9N2qemuALJOfiapejej/I9UA/kdE3lOs287y/czu2276cyMKRlQ1dAMRlZCI/BnRu5e8p93l3wUwRVUPGLQRkbkALgfwTVUtdAKWDPBnQkSF4iN9RGUkPv+oAcDiDJ8ej8yP8gHRxPFOFH+AiDqPPxMiKgg3fUTlZSSiyd82m7t44reh/eXx5wYgOmfzLu3c+xFTkfFnQkSdwffeJSovqXeNab+5q0M0LX7Apk9VNyP76y9SAPyZEFFn8Jw+IiIiojLAp3eJiIiIygA3fURERERlgJs+IiIiojLATR8RERFRGeCmj4iIiKgMcNNHREREVAa46SMiIiIqA9z0EREREZUBbvqIiIiIygA3fURERERlgJs+IiIiojLATR8RERFRGeCmj4iIiKgMcNNHREREVAa46SMiIiIqA9z0EREREZUBbvqIiIiIygA3fURERERlgJs+IiIiojLATR8RERFRGeCmj4iIiKgMcNNHREREVAa46SMiIiIqA9z0EREREZUBbvqIiIiIygA3fURERERlgJs+IiIiojLATR8RERFRGagKHVCo/v37a21tbegMIiKikli4cOFrqjogdAf5527TV1tbiwULFoTOICIiKgkRWRO6gZKBT+8SERERlQFu+oiIiIjKADd9RERERGWAmz4iIiKiMsBNHxEREVEZ4KaPiIiIqAxw00dERERUBrjpIyIiIioD3PQRERERlQHTTZ+ITBGRFSKySkRmZ/j8RSKyWUSej9elXfl+TXfORO3VVai4UVB7dRWa7pzZlZuz1dQE1NYCFRXRr01NoYuy89LqpRNgqxW22mCrCVfHLEoGVTVZACoBvARgJIBqAIsB1Le7zkUA/l8htztx4kTN5KE7GrXmOihu3L9qroM+dEdjxusH9dBDqjU1qsD+VVMTXd7deGn10qnKVitstcFWE4UcswAsUKNjNVd5LVFVk82kiEwGcKOqvjf++Jp4k/nVtOtcBGCSql6Z7+1OmjRJM733bu3VVVjTZ98Blw/fXonVt7UU3G+qthZYk+GtFIcPB1avLnVNbl5avXQCbLXCVhtsNVHIMUtEFqrqpFK1UXJZPr07GMC6tI9fiS9r7xwRWSIij4rI0Ew3JCIzRGSBiCzYvHlzxm+29qAD/+PJdXlQa9cWdnlIXlq9dAJstcJWG2w14eqYRYkRepDjcQC1qnosgF8BuD/TlVT1blWdpKqTBgwYkPGGhr1dWdDlQQ0bVtjlIXlp9dIJsNUKW22w1YSrYxYlhuWmbz2A9EfuhsSX/ZOqblHV3fGH9wCY2NlvdsvIGajZ2/aymr3R5d3OLbcANTVtL6upiS7vbry0eukE2GqFrTbYasLVMYuSw+pkQQBVAF4GMAL7Bzka2l1nYNrvPwTguY5uN9sgh6pq4+frtfKLUNwArfwitPHz9VmvG1xjo2plpSoQ/drYDQdOUry0eulUZasVttpgq4l8j1ngIAdXkZbtjQPvB7AS0RTvdfFlNwGYGv/+qwCWxhvC3wIY09Ftcnq3xLy0eulUZasVttpgqwlO73KFWGbTu1Y4vVtiXlq9dAJstcJWG2w1weldCiExm76KGwUqB15fFGi9sZv9GSsqon+DticCtLaWvicXL61eOgG2WmGrDbaaKOSYxU0fFUvo6d2icTUJ5WjCzE2rl06ArVbYaoOtJlwdsygxErPpczUJ5WjCzE2rl06ArVbYaoOtJlwdsyg5Qp9UWOji9G4AXlq9dKqy1QpbbbDVBKd3uUq9ggcUuji9W2JeWr10qrLVClttsNUEp3e5QqzEDHJweteIl1YvnQBbrbDVBltNcHqXQkjMpo/Tu0a8tHrpBNhqha022GqC07sUQmIGOVxNQjmaMHPT6qUTYKsVttpgqwlXxyxKjMRs+lxNQjmaMHPT6qUTYKsVttpgqwlXxyxKjtAnFRa6OL0bgJdWL52qbLXCVhtsNcHpXa5Sr+ABhS5O75aYl1YvnapstcJWG2w1weldrhArMYMcnN414qXVSyfAVitstcFWE5zepRASs+nj9K4RL61eOgG2WmGrDbaa4PQuhZCYQQ5Xk1COJszctHrpBNhqha022GrC1TGLEiMxmz5Xk1COJszctHrpBNhqha022GrC1TGLkiP0SYWFLk7vBuCl1UunKlutsNUGW01weper1Ct4QKGL07sl5qXVS6cqW62w1QZbTXB6lyvESswgB6d3jXhp9dIJsNUKW22w1QSndymExGz6OL1rxEurl06ArVbYaoOtJji9SyEkZpDD1SSUowkzN61eOgG2WmGrDbaacHXMosRIzKbP1SSUowkzN61eOgG2WmGrDbaacHXMouQIfVJhoYvTuwF4afXSqcpWK2y1wVYTnN7lKvUKHlDo4vRuiXlp9dKpylYrbLXBVhOc3uUKsRIzyMHpXSNeWr10Amy1wlYbbDXB6V0KITGbPk7vGvHS6qUTYKsVttpgqwlO71IIiRnkcDUJ5WjCzE2rl06ArVbYaoOtJlwdsygxErPpczUJ5WjCzE2rl06ArVbYaoOtJlwdsyg5Qp9UWOji9G4AXlq9dKqy1QpbbbDVBKd3uUq9ggcUuji9W2JeWr10qrLVClttsNUEp3e5QqzEDHJweteIl1YvnQBbrbDVBltNcHqXQkjMpo/Tu0a8tHrpBNhqha022GqC07sUQmIGOVxNQjmaMHPT6qUTYKsVttpgqwlXxyxKjMRs+lxNQjmaMHPT6qUTYKsVttpgqwlXxyxKjtAnFRa6OL0bgJdWL52qbLXCVhtsNcHpXa5Sr+ABhS5O75aYl1YvnapstcJWG2w1weldrhArMYMcnN414qXVSyfAVitstcFWE5zepRASs+nj9K4RL61eOgG2WmGrDbaa4PQuhZCYQQ5Xk1COJszctHrpBNhqha022GrC1TGLEiMxmz5Xk1COJszctHrpBNhqha022GrC1TGLkiP0SYWFLk7vBuCl1UunKlutsNUGW01weper1Ct4QKGL07sl5qXVS6cqW62w1QZbTXB6lyvESswgB6d3jXhp9dIJsNUKW22w1QSndymExGz6OL1rxEurl06ArVbYaoOtJji9SyEkZpDD1SSUowkzN61eOgG2WmGrDbaacHXMosRIzKbP1SSUowkzN61eOgG2WmGrDbaacHXMouQIfVJhoYvTuwF4afXSqcpWK2y1wVYTnN7lKvUKHlDo4vRuiXlp9dKpylYrbLXBVhOc3uUKsRIzyMHpXSNeWr10Amy1wlYbbDXB6V0KITGbPk7vGvHS6qUTYKsVttpgqwlO71IIiRnkcDUJ5WjCzE2rl06ArVbYaoOtJlwdsygxTDd9IjJFRFaIyCoRmZ3jeueIiIpIp/8l42oSytGEmZtWL50AW62w1QZbTbg6ZlFyWJ0sCKASwEsARgKoBrAYQH2G6/UF8DsAzwGY1NHtcno3AC+tXjpV2WqFrTbYaoLTu1ylXpaP9J0AYJWqvqyqewD8EMDZGa53M4CvAdjVlW/WdOdM3N+jGfsqAQiwrxK4v0czmu6c2ZWbtdHUBNx/P7AvPol3377o46amsF2ZeGn10gmw1QpbbbDVhKtjFiWG2SCHiHwEwBRVvTT++AIAJ6rqlWnXOR7Adap6jog8BeBzqnrglEYaTu+WmJdWL50AW62w1QZbTXB6l0IINsghIhUA/gvAf+Rx3RkiskBEFmzevDnjddYedOB/PLkuD2rt2sIuD8lLq5dOgK1W2GqDrSZcHbMoMSw3fesBDE37eEh8WUpfAMcAeEpEVgM4CcD8TMMcqnq3qk5S1UkDBgzI+M1cTUI5mjBz0+qlE2CrFbbaYKsJV8csSgzLTd9fAIwWkREiUg1gGoD5qU+q6jZV7a+qtapai2iQY2pHT+9m42oSytGEmZtWL50AW62w1QZbTbg6ZlFyWE6JAHg/gJWIpniviy+7CdHmrv11nwKnd7snL61eOlXZaoWtNthqgtO7XKVewQMKXXzv3RLz0uqlU5WtVthqg60m+N67XCFWYt6GjdO7Rry0eukE2GqFrTbYaoLTuxRCYjZ9fO9dI15avXQCbLXCVhtsNcH33qUQ+N67ITiaMHPT6qUTYKsVttpgqwlXxyxKjMRs+lxNQjmaMHPT6qUTYKsVttpgqwlXxyxKjtAnFRa6OL0bgJdWL52qbLXCVhtsNcHpXa5Sr+ABhS5O75aYl1YvnapstcJWG2w1weldrhArMYMcnN414qXVSyfAVitstcFWE5zepRASs+nj9K4RL61eOgG2WmGrDbaa4PQuhZCYQQ5Xk1COJszctHrpBNhqha022GrC1TGLEiMxmz5Xk1COJszctHrpBNhqha022GrC1TGLkiP0SYWFLk7vBuCl1UunKlutsNUGW01weper1Ct4QKGL07sl5qXVS6cqW62w1QZbTXB6lyvESswgB6d3jXhp9dIJsNUKW22w1QSndymExGz6OL1rxEurl06ArVbYaoOtJji9SyEkZpDD1SSUowkzN61eOgG2WmGrDbaacHXMosRIzKbvlpEzUN3uWdzqlm46CXXLLUB1ddvLqqu75YSZm1YvnQBbrbDVBltNuDpmUWIkZtMHAO0f1O9mT+q21f4piO78NLuXVi+dAFutsNUGW024OmZRIiTmnD4Ochjx0uqlE2CrFbbaYKsJDnJQCInZ9HGQw4iXVi+dAFutsNUGW01wkINCSMzTu65OinV0srGbVi+dAFutsNUGW024OmZRYiRm0+fqLW0cvVWQm1YvnQBbrbDVBltNuDpmUXKEfnXoQhffhi0AL61eOlXZaoWtNthqgm/DxlXqFTyg0MW3YSsxL61eOlXZaoWtNthqgm/DxhViJWaQg9O7Rry0eukE2GqFrTbYaoLTuxRCYjZ9nN414qXVSyfAVitstcFWE5zepRASM8jhahLK0YSZm1YvnQBbrbDVBltNuDpmUWIkZtPnahLK0YSZm1YvnQBbrbDVBltNuDpmUXKEPqmw0MXp3QC8tHrpVGWrFbbaYKsJTu9ylXoFDyh0cXq3xLy0eulUZasVttpgqwlO73KFWIkZ5OD0rhEvrV46AbZaYasNtprg9C6FkJhNH6d3jXhp9dIJsNUKW22w1QSndymExAxyuJqEcjRh5qbVSyfAVitstcFWE66OWZQYidn0uZqEcjRh5qbVSyfAVitstcFWE66OWZQcoU8qLHRxejcAL61eOlXZaoWtNthqgtO7XKVewQMKXZzeLTEvrV46Vdlqha022GqC07tcIVZiBjk4vWvES6uXToCtVthqg60mOL1LISRm08fpXSNeWr10Amy1wlYbbDXB6V0KITGDHK4moRxNmLlp9dIJsNUKW22w1YSrYxYlRmI2fa4moRxNmLlp9dIJsNUKW22w1YSrYxYlR+iTCgtdnN4NwEurl05Vtlphqw22muD0LlepV/CAQhend0vMS6uXTlW2WmGrDbaa4PQuV4iVmEEOTu8a8dLqpRNgqxW22mCrCU7vUgiJ2fRxeteIl1YvnQBbrbDVBltNcHqXQkjMIIerSShHE2ZuWr10Amy1wlYbbDXh6phFiZGYTZ+rSShHE2ZuWr10Amy1wlYbbDXh6phFyRH6pMJCF6d3A/DS6qVTla1W2GqDrSY4vctV6hU8oNDF6d0S89LqpVOVrVbYaoOtJji9yxViJWaQg9O7Rry0eukE2GqFrTbYaoLTuxRCYjZ9nN414qXVSyfAVitstcFWE5zepRASM8jhahLK0YSZm1YvnQBbrbDVBltNuDpmUWIkZtPnahLK0YSZm1YvnQBbrbDVBltNuDpmUXJYnjAIYAqAFQBWAZid4fNXAHgBwPMAngFQ39Ftcno3AC+tXjpV2WqFrTbYaoLTu1ylXnY3DFQCeAnASADVABa333yVS8QAACAASURBVNQBODjt91MB/Lyj2+X0bol5afXSqcpWK2y1wVYTnN7lCrHMBjlEZDKAG1X1vfHH18SPLH41y/U/CuBCVX1frtvl9G6JeWn10gmw1QpbbbDVBKd3KYQqw9seDGBd2sevADix/ZVE5N8BfBbRo4H/mumGRGQGgBkAMCzLCblrDzrwP55clwe1dm1hl4fkpdVLJ8BWK2y1wVYTro5ZlBjBBzlU9TuqehSA/wTwhSzXuVtVJ6nqpAEDBmS8HVeTUI4mzNy0eukE2GqFrTbYasLVMYsSw3LTtx7A0LSPh8SXZfNDAB/s7DdzNQnlaMLMTauXToCtVthqg60mXB2zKDmsThZE9NTxywBGYP8gR0O764xO+/1ZyONkVU7vBuCl1UunKlutsNUGW01weper1Mv2xoH3A1iJaIr3uviymwBMjX//3wCWInrJlt+23xRmWpzeLTEvrV46Vdlqha022GqC07tcIVZi3oaN07tGvLR66QTYaoWtNthqgtO7FEJiNn18710jXlq9dAJstcJWG2w1wffepRDyGuQQkXeIyIdE5N9F5BIROUFEgk/+pnM1CeVowsxNq5dOgK1W2GqDrSZcHbMoMXJu3ETkdBH5BYCfAngfgIEA6hG9tMoLIvIlETnYPrNjriahHE2YuWn10gmw1QpbbbDVhKtjFiVHrhP+ANwGYFiWz1UheomVc0p5EiKndwPw0uqlU5WtVthqg60mOL3LVeoVPKDQxendEvPS6qVTla1W2GqDrSY4vcsVYuU1yCEinwFwH4C3ANwDYAKA2ar6S5OHH3Pg9G6JeWn10gmw1QpbbbDVBKd3KYR8N32LVXW8iLwXwOUAvgjgQVU93jqwPU7vlpiXVi+dAFutsNUGW01wepdCqMrzeqm/mu9HtNlbKiIZ/rqGM+ztyoz/auqWk1DDhmX+12g3nDBz0+qlE2CrFbbaYKuJrh6zFi5c+I6qqqp7ABwD27dUJT9aAbzY0tJy6cSJE1/NdIV8N30LReSXiN5S7RoR6RvfeLdxy8gZmLH+Tuzosf+ybjsJdcstwIwZwI4d+y/rphNmblq9dAJstcJWG2w10dVjVlVV1T1HHnnk2AEDBrxeUVHRzZ7OohBaW1tl8+bN9Zs2bboHwNRM18n3XwefBDAbwDtVdQei99K9uDiZxTG98Q58Ym89KvcBUKByH/CJvfWY3nhH6LQDTZ8OfOITQGX8L7rKyujj6dPDdmXipdVLJ8BWK2y1wVYTRThmHTNgwIA3ueGjlIqKCh0wYMA2RI/+Zr5OPjekqq0AWgCcIiIfBnAqgFFFqSySpjtn4v4ezdhXCUCAfZXA/T2a0XTnzNBpB2pqAu6/H9gXP7S/b1/0cVNT2K5MvLR66QTYaoWtNthqogjHrApu+Ki9+O9E1r1dvoMc9wI4FsBS7H9aV1X1kmJEFoLTuyXmpdVLJ8BWK2y1wVYTXZ3eXbx48erx48e/ZltJHi1evLj/+PHjazN9Lt+nd09S1Umq+glVvTheJd/w5bL2oAP/48l1eVBr1xZ2eUheWr10Amy1wlYbbDXh6piVRU1NzYSOrnP++ecPX7hwYS8AmD179pHpn5swYcKYzn6PtWvXVn3gAx8YOXTo0GMaGhrGnnrqqaOWLFnSM9/2bJ544om+p59++igAaGpq6nfttdceCQAPPvjgIak/BwBcddVVgx577LG+Xf1+pZbvpu9ZEak3LekiV+9j6Oj9Id20eukE2GqFrTbYaqLUx6y5c3HYoEEYV1GBiYMGYdzcuTjM5Bu18/DDD6+ZOHHiLgCYM2fOwPTPLVq0aHlnbrO1tRVTp04ddcopp7y1bt26F5cuXbrs1ltvXb9hw4YeHX91/qZPn77tK1/5yiYAeOyxxw5ZsmRJ79TnvvWtb2344Ac/+FYxv18p5LvpewDRxm+FiCwRkRdEZIllWKFcvY+ho/eHdNPqpRNgqxW22mCriVIes+bOxWGzZmH4xo2oVgU2bkT1rFkYXqyN3xNPPNH3hBNOOHrKlCkjR4wY0TB16tQRrfHrIp5wwglH/+53v6uZOXPm4N27d1eMGTOmfurUqSOA/Y/ibdu2rWLy5Ml19fX1Y+vq6uofeuihQzr6flVVVfr5z39+c+qyyZMn75wyZcr21tZWXH755UNGjx7dUFdXVz9v3rxDO2p89NFHDx4xYkRDfX392EcfffSf33vOnDmHX3jhhcN+9atfHfTrX//6kC984QtDxowZU7906dKe55xzTu199913KAD85Cc/6Tt27Nj6urq6+nPPPbd2586dAgCDBw8eN2vWrEGpP9eiRYt6IbB8N33fBXABgCkAzgLwgfjXboPTu0a8tHrpBNhqha022GqilMesm27C4F272h7vd+1CxU03YXCxvseyZct6f+c731m3atWqpWvXru35q1/9qk/65++44471PXv2bF2+fHnz/Pnz/57+uZqamtaf/vSnq5qbm5c9/fTTK6+99tohrTleTHvJkiW9x48fvyPT5x544IFDXnjhhd7Lli1b+uSTT668/vrrh6xZs6ZHtsYdO3bIlVdeWTt//vxVL7744rJXX331gEcL3/3ud7995plnvvHlL3/5leXLlzc3NDTsTn1ux44dcvnll494+OGHX1q5cmVzS0sLbrvttgGpz/fv37+lubl52SWXXLL51ltvPSLvO9RIvpu+zao6X1X/rqprUsu0rECc3jXipdVLJ8BWK2y1wVYTpTxmbdqE6kIu74xx48a9fdRRR+2trKxEQ0PDjpdeeinv225tbZWrrrpqSF1dXf3pp59e9+qrr1a/8sor+b6OcBu///3v+5533nlbq6qqMHTo0JYTTzxx+zPPPFOTrfH555/vNWTIkN3jxo3bXVFRgenTp28p5PstXry415AhQ3Yfe+yxuwHgoosu2vLMM8/881y/j33sY68DwAknnLBj3bp1XT7nsKvy3fQtEpHvi8hHReTDqWVaVqDrXr67zYtcAsCOHtHl3c5117V98VAg+vi668L05OKl1UsnwFYrbLXBVhOlPGYdeST2FHJ5Z/Ts2fOfLwVSWVmJlpaWvN+166677jpsy5YtVS+88MKy5cuXNx9++OF7d+7cmXV/Mm7cuJ2LFy+uyfZ5i8bO6tWrlwJAVVWVluL7dSTfTV9vALsBvAfR07qpp3i7DVeTUI4mzNy0eukE2GqFrTbYaqKUx6zrr8f6Xr3avotWr15ovf56rC/6N8uhqqpKd+/efcDGZ9u2bZX9+/ff27NnT3388cf7btiwIeejhGedddZbe/bskW984xv9U5f96U9/6v3zn/+8zymnnPLWo48+elhLSws2bNhQ9ec//7nPySef/Ha22zruuON2rV+/vnrp0qU9AeCHP/xhxvMc+/Tps+/NN988YM80fvz4XevXr69+8cUXewLAAw88cPjJJ5/cbQc88n1x5oszrG71ki2c3jXipdVLJ8BWK2y1wVYTpTxmXXEFtt5+O9YMHIg9IsDAgdhz++1Yc8UV2Fr0b5bD9OnTN48dO/afgxwpl1566dbFixcfVFdXV3///fcfPmLEiF25bqeiogLz589/6Te/+c3BQ4cOPWbUqFEN//mf/zl48ODBey+44II3Ghoado4dO7bhtNNOq/vSl770yrBhw7K+WG9NTY1++9vfXvOBD3xgVH19/dj+/ftnvO706dO3zpkz58ixY8fWpzaIqa+fO3fu6nPPPfeourq6+oqKCnzuc5/bnOk2uoOcL84sIl8AcIeqZvyLISL/CqBGVZ8w6jtAthdnbrpzZsb3Mbx7cGP3G+Zoasr8/pB33939Tjj20uqlE2CrFbbaYKuJQo5ZfHFmKkRXXpz5BQCPi8iTInKbiHxeRK4XkQdF5AVET/P+qci9ncLpXSNeWr10Amy1wlYbbDXh6phFiZFz06eqP1HVdwG4AtFbsFUCeBPAQwBOUNVZqtotHsbk9K4RL61eOgG2WmGrDbaacHXMosTI6713uxO+926JeWn10gmw1QpbbbDVBN97l6wU4713uz1O7xrx0uqlE2CrFbbaYKsJV8csSozEbPo4vWvES6uXToCtVthqg60mXB2zKDESs+nje+8a8dLqpRNgqxW22mCrCVfHLEqMvDZ9IlIXT/C+GH98bPxyLt2Gq0koRxNmblq9dAJstcJWG2w14eqYlUVlZeXEMWPG1B999NH19fX1Y3/1q18dZP09Bw8ePG7jxo2deou2XHbv3i0zZ84cPHz48GPq6+vHHnfccWMeeeSRg4tx2zU1NRMAYPXq1T2mTJkyEgD++Mc/9n744Yf7pa7T1NTU79prrz2yGN8vl3wf6ZsH4BoAewFAVZcAmGYV1RmuJqEcTZi5afXSCbDVClttsNVEqY9Zc/8y97BB3xw0ruJLFRMHfXPQuLl/mZvxnScK0bNnz9bly5c3r1ixovnmm29ef+211w4pRmsxtbTkN8g5a9asQZs2beqxfPnypc3Nzcsef/zxVW+++WZRn2uvra3d+/Of//xlAFiwYEHNT3/6039u+qZPn77tK1/5yqZifr9M8t301ajqn9td1q1GYvneu0a8tHrpBNhqha022GqilMesuX+Ze9isX84avnH7xmqFYuP2jdWzfjlreDE2finbtm2r7NevXwsAtLa24vLLLx8yevTohrq6uvp58+YdCgBPPPFE39NPP31U6msuvPDCYXPmzDkciB7BmzVr1qD6+vqxdXV19YsWLeoFAJs2bap817veNXrUqFEN559//vD0Vxw588wzj2poaBg7atSohvS3ZKupqZlw2WWXDTn66KPrZ8+ePfDMM888KvW5H//4xwe/+93v/ufHAPDWW29VfP/73x9wzz33rO3du7cCwNChQ1suvfTS14HovYHr6urqR48e3dDY2Dg4/ft86lOfGnz00UfXjx8/fsy6deuqAGD58uXVxx133Ji6urr6T3/604NS11+xYkX16NGjG3bt2iVf/epXBz3++OOHjhkzpn7evHmHzpkz5/ALL7xwWOp6J510Ul1dXV395MmT6/72t79VA8A555xTe9FFFw2dMGHCmCFDhoy77777Di3055Tvpu81ETkKgAKAiHwEwMZCv5klV5NQjibM3LR66QTYaoWtNthqopTHrJt+d9PgXS272hzvd7XsqrjpdzcNzvY1+di9e3fFmDFj6keMGNHwmc98ZvgNN9ywEQAeeOCBQ1544YXey5YtW/rkk0+uvP7664esWbOmR0e3179//5bm5uZll1xyyeZbb731CACYPXv2oMmTJ29ftWrV0g996ENvbNy48Z/vy9vU1LR66dKly55//vnmu+6664hNmzZVAsDOnTsrTjzxxLdXrFjR/PWvf33jSy+91GvDhg1VAHDvvfcefvHFF7d5qZvm5uaeAwcO3HPYYYe1eX9iIHpK9sYbbxz81FNPrWxubl66aNGigx588MFDUt9n8uTJ21esWNE8efLk7d/+9rcHAMDMmTOHXXrppZtXrlzZPHDgwL3tb7NXr156zTXXbDjrrLNeX758efNll132evrnGxsbh02fPn3LypUrm88///wtjY2NQ1Of+8c//tFjwYIFy3/yk5/87YYbbij455fvpu/fAdwFYIyIrAdwFYDGQr+ZJVeTUI4mzNy0eukE2GqFrTbYaqKUx6xN2zdVF3J5vlJP7/79739f+uMf//hvF1988YjW1lb8/ve/73veeedtraqqwtChQ1tOPPHE7c8880xNR7f3sY997HUAOOGEE3asW7euJwA899xzfS+55JItADBt2rRtBx988D93xV/72teOOProo+snTpw4dtOmTT2WLl3aCwAqKytx0UUXvQ5E79N73nnnbZk3b95hr732WuVf//rXPueee+62fP+MzzzzzEEnnXTSW4MGDWrp0aMHzj///K1PP/10HwDo0aOHTps2bRsATJw48e01a9ZUA8Bf//rXPpdddtlWALj88su35Pu9UhYtWnTQjBkztgJAY2Pj1oULF/ZJfW7q1KlvVFZWYuLEibu2bNnS4Ua6vbw2far6sqqeCWAAgDGq+v+p6upCv5klV5NQjibM3LR66QTYaoWtNthqopTHrCP7HLmnkMs748wzz3z79ddfr8o1ZNGjRw9tbd3/YNru3bsl/fO9evVSAKiqqtKWlhZBDk888UTfp59+uu+CBQuWr1ixonns2LE7d+7cWQEA1dXVrVVV+zMaGxu3PPLII4d/97vfPeyss856vUePtnul+vr63Rs3bqzeunVrQa9oUlVVpRUVFanfI725oqLC5J0vUvcRAHTmzTXynd49REQ+DeBmALeIyBwRmVPwdzPkahLK0YSZm1YvnQBbrbDVBltNlPKYdf0p16/vVdWrzVOXvap6tV5/yvXri/U9Fi1a1Ku1tRVHHHFEyymnnPLWo48+elhLSws2bNhQ9ec//7nPySef/PZRRx21e9WqVb137twpr732WuUzzzzT4XTsSSed9Nb3vve9wwHgkUceOTg1XPHGG29U9uvXb1/fvn1bFy1a1Gvx4sVZJ4dra2v3HnHEEXu/+c1vDpwxY8YB72LSt2/f1mnTpr02Y8aMYbt27RIA2LBhQ9W999576Mknn/z2n/70p74bN26samlpwY9+9KPDTjvttO25mo8//vjt8+bNOwwA5s2bd3im6xx88MH7tm/fnnEPNmHChLfvueeeQ4HofMJJkybl/H6FyHdX+zMAtQBeALAwbXUbnN414qXVSyfAVitstcFWE6U8Zl3xziu23v6e29cM7DNwj0AwsM/APbe/5/Y1V7zziq1dud3UOX1jxoypnzZt2sg777xzdVVVFS644II3Ghoado4dO7bhtNNOq/vSl770yrBhw1pGjRq196yzznp9zJgxDWefffbIhoaGHR19j1tvvXXDH/7whz6jRo1q+N///d9DBw4cuAcAzjnnnG0tLS0ycuTIhquvvnrw+PHj3851O9OmTdsycODAPccff/yuTJ//1re+tb5///4tdXV1DaNHj26YMmXKqH79+u0bPnz43htuuGH9qaeeWjd27NiG8ePHv/3xj3/8jVzf64477lh79913v6Ourq5+/fr1GZ+Cfd/73vfWypUre6cGOdI/N3fu3LUPPvhg/7q6uvof/OAHh99xxx3rct9L+cvrvXdF5K+qenyxvmlX8L13S8xLq5dOgK1W2GqDrSb43ruldeGFFw6bMGHCjlmzZiX+Psv13rv5vsDhgyJyGYAnAOxOXaiqXfpXQjFxeteIl1YvnQBbrbDVBltNuDpmOdfQ0DC2d+/erXfddVfRHjHzKt+nd/cAuA3As9j/1O6BD7cFxOldI15avXQCbLXCVhtsNeHqmOXc0qVLly1YsGBF6jX4ylm+m77/ADBKVWtVdUS8RlqGFYrTu0a8tHrpBNhqha022GqiCMes1tbW1pwTrlR+4r8TB7zeYEq+m75VADo84TIkTu8a8dLqpRNgqxW22mCriSIcs17cvHlzP278KKW1tVU2b97cD8CL2a6T7yDHjwE0APgt2p7T9+kidBYk2yBH050zMWP9nW3e1qZmL3D34Mbut/FragJmzGj7dkE1NcDdd3e//zl5afXSCbDVClttsNVEIcesTIMcCxcufEdVVdU9AI5B/g/gULK1AnixpaXl0okTJ76a6Qr5bvo+kelyVb2/a32F4/RuiXlp9dIJsNUKW22w1URXp3eJOiOvTV93km3TV3GjQDM8yC0KtN7Yzf6MFRVApvtdBGjN+lR8GF5avXQCbLXCVhtsNVHIMYubPiqWnA8Ji8gj8a8viMiS9qs0iflxNQnlaMLMTauXToCtVthqg60mXB2zKDE6Og/gM/GvHwBwVobVbdwycgaq2z2LW93Sjad3q9u9z3V1dbecMHPT6qUTYKsVttpgqwlXxyxKjJybPlXdGP92pqquSV8Aut37m7V/UL+bPanbVvunILrz0+xeWr10Amy1wlYbbDXh6phFidDpt2ETkSWqeqxZWRYc5CgxL61eOgG2WmGrDbaa4CAHhZBz0ycijYge0RsJ4KW0T/UF8AdV/bht3oE4yFFiXlq9dAJstcJWG2w1wUEOCqGjc/q+j+jcvfloey7fxBAbvlxcnRTr6GRjN61eOgG2WmGrDbaacHXMosTo6Jy+baq6WlU/2u6cvq2lCswX34bNiJdWL50AW62w1QZbTbg6ZlFyqKrZAjAFwApEb+M2O8PnPwugGcASAE8CGN7RbU6cOFGzafx8vVZ+EYoboJVfhDZ+vj7rdYNrbFStrFQFol8bG0MXZeel1UunKlutsNUGW03ke8wCsEANj9Vc5bPsbhioRHQe4EgA1QAWA6hvd53TAdTEv28E8HBHt5tt0/fQHY1acx0UN+5fNddBH7qjG/4H/9BDqjU10d2fWjU10eXdjZdWL52qbLXCVhtsNVHIMYubPq5iLbN35BCRyQBuVNX3xh9fEz+y+NUs158A4P+p6rty3S6nd0vMS6uXToCtVthqg60mOL1LIVi+SfNgAOvSPn4lviybTwL4v0yfEJEZIrJARBZs3rw54xevPejA/3hyXR7U2rWFXR6Sl1YvnQBbrbDVBltNuDpmUWJYbvryJiIfBzAJwG2ZPq+qd6vqJFWdNGDAgIy34WoSytGEmZtWL50AW62w1QZbTbg6ZlFiWG761gMYmvbxkPiyNkTkTADXAZiqqrs7+81cTUI5mjBz0+qlE2CrFbbaYKsJV8csSg6rkwUBVAF4GcAI7B/kaGh3nQmIhj1G53u7nN4NwEurl05Vtlphqw22muD0Llepl+2NA+8HsDLe2F0XX3YTokf1AODXAP4B4Pl4ze/oNjm9W2JeWr10qrLVClttsNUEp3e5Qiyz6V0rnN4tMS+tXjoBtlphqw22muD0LoWQmE0f33vXiJdWL50AW62w1QZbTfC9dymEbjG9WwyuJqEcTZi5afXSCbDVClttsNWEq2MWJUZiNn2uJqEcTZi5afXSCbDVClttsNWEq2MWJUfokwoLXZzeDcBLq5dOVbZaYasNtprg9C5XqVfwgEIXp3dLzEurl05Vtlphqw22muD0LleIlZhBDk7vGvHS6qUTYKsVttpgqwlO71IIidn0cXrXiJdWL50AW62w1QZbTXB6l0JIzCCHq0koRxNmblq9dAJstcJWG2w14eqYRYmRmE2fq0koRxNmblq9dAJstcJWG2w14eqYRckR+qTCQhendwPw0uqlU5WtVthqg60mOL3LVeoVPKDQxendEvPS6qVTla1W2GqDrSY4vcsVYiVmkIPTu0a8tHrpBNhqha022GqC07sUQmI2fZzeNeKl1UsnwFYrbLXBVhOc3qUQEjPI4WoSytGEmZtWL50AW62w1QZbTbg6ZlFiJGbT52oSytGEmZtWL50AW62w1QZbTbg6ZlFyhD6psNDF6d0AvLR66VRlqxW22mCrCU7vcpV6BQ8odHF6t8S8tHrpVGWrFbbaYKsJTu9yhViJGeTg9K4RL61eOgG2WmGrDbaa4PQuhZCYTR+nd414afXSCbDVClttsNUEp3cphMQMcriahHI0Yeam1UsnwFYrbLXBVhOujlmUGInZ9LmahHI0Yeam1UsnwFYrbLXBVhOujlmUHKFPKix0cXo3AC+tXjpV2WqFrTbYaoLTu1ylXsEDCl2c3i0xL61eOlXZaoWtNthqgtO7XCFWYgY5OL1rxEurl06ArVbYaoOtJji9SyEkZtPH6V0jXlq9dAJstcJWG2w1weldCiExgxyuJqEcTZi5afXSCbDVClttsNWEq2MWJUZiNn2uJqEcTZi5afXSCbDVClttsNWEq2MWJUfokwoLXZzeDcBLq5dOVbZaYasNtprg9C5XqVfwgEIXp3dLzEurl05Vtlphqw22muD0LleIlZhBDk7vGvHS6qUTYKsVttpgqwlO71IIidn0cXrXiJdWL50AW62w1QZbTXB6l0JIzCCHq0koRxNmblq9dAJstcJWG2w14eqYRYmRmE2fq0koRxNmblq9dAJstcJWG2w14eqYRckR+qTCQhendwPw0uqlU5WtVthqg60mOL3LVeoVPKDQxendEvPS6qVTla1W2GqDrSY4vcsVYiVmkIPTu0a8tHrpBNhqha022GqC07sUQmI2fZzeNeKl1UsnwFYrbLXBVhOc3qUQEjPI4WoSytGEmZtWL50AW62w1QZbTbg6ZlFiJGbT52oSytGEmZtWL50AW62w1QZbTbg6ZlFyhD6psNDF6d0AvLR66VRlqxW22mCrCU7vcpV6BQ8odHF6t8S8tHrpVGWrFbbaYKsJTu9yhViJGeTg9K4RL61eOgG2WmGrDbaa4PQuhZCYTR+nd414afXSCbDVClttsNUEp3cphMQMcriahHI0Yeam1UsnwFYrbLXBVhOujlmUGInZ9LmahHI0Yeam1UsnwFYrbLXBVhOujlmUHKFPKix0cXo3AC+tXjpV2WqFrTbYaoLTu1ylXsEDCl2c3i0xL61eOlXZaoWtNthqgtO7XCFWYgY5OL1rxEurl06ArVbYaoOtJji9SyEkZtPH6V0jXlq9dAJstcJWG2w1weldCiExgxyuJqEcTZi5afXSCbDVClttsNWEq2MWJYbppk9EpojIChFZJSKzM3z+FBH5q4i0iMhHuvK9XE1COZowc9PqpRNgqxW22mCrCVfHLEoOq5MFAVQCeAnASADVABYDqG93nVoAxwJ4AMBH8rldTu8G4KXVS6cqW62w1QZbTXB6l6vUy/KRvhMArFLVl1V1D4AfAji73YZztaouAdDlky2a7pyJ+3s0Y18lAAH2VQL392hG050zu3rTxdfUBNx/P7AvPol3377o46amsF2ZeGn10gmw1QpbbbDVhKtjFiWG2SBH/HTtFFW9NP74AgAnquqVGa77PQBPqOqjHd0up3dLzEurl06ArVbYaoOtJji9SyG4GOQQkRkiwQRC/QAADjxJREFUskBEFmzevDnjddYedOB/PLkuD2rt2sIuD8lLq5dOgK1W2GqDrSZcHbMoMSw3fesBDE37eEh8WcFU9W5VnaSqkwYMGJDxOq4moRxNmLlp9dIJsNUKW22w1YSrYxYlhuWm7y8ARovICBGpBjANwHyrb+ZqEsrRhJmbVi+dAFutsNUGW024OmZRclhOiQB4P4CViKZ4r4svuwnA1Pj37wTwCoC3AWwBsLSj2+T0bgBeWr10qrLVClttsNUEp3e5Sr2CBxS6+N67Jeal1UunKlutsNUGW03wvXe5QqzEvA0bp3eNeGn10gmw1QpbbbDVBKd3KYTEbPr43rtGvLR66QTYaoWtNthqgu+9SyG4eMmWfLiahHI0Yeam1UsnwFYrbLXBVhOujlmUGInZ9LmahHI0Yeam1UsnwFYrbLXBVhOujlmUHKFPKix0cXo3AC+tXjpV2WqFrTbYaoLTu1ylXsEDCl2c3i0xL61eOlXZaoWtNthqgtO7XCFWYgY5OL1rxEurl06ArVbYaoOtJji9SyEkZtPH6V0jXlq9dAJstcJWG2w1weldCiExgxyuJqEcTZi5afXSCbDVClttsNWEq2MWJUZiNn23jJyB6nbP4la3dNNJqFtuAaqr215WXd0tJ8zctHrpBNhqha022GrC1TGLEiMxmz4AaP+gfjd7Uret9k9BdOen2b20eukE2GqFrTbYasLVMYsSITHn9HGQw4iXVi+dAFutsNUGW01wkINCSMymj4McRry0eukE2GqFrTbYaoKDHBRCYp7edXVSrKOTjd20eukE2GqFrTbYasLVMYsSIzGbPldvaePorYLctHrpBNhqha022GrC1TGLkiP0q0MXuvg2bAF4afXSqcpWK2y1wVYTfBs2rlKv4AGFLr4NW4l5afXSqcpWK2y1wVYTfBs2rhArMYMcnN414qXVSyfAVitstcFWE5zepRASs+nj9K4RL61eOgG2WmGrDbaa4PQuhZCYQQ5Xk1COJszctHrpBNhqha022GrC1TGLEiMxmz5Xk1COJszctHrpBNhqha022GrC1TGLkiP0SYWFLk7vBuCl1UunKlutsNUGW01weper1Ct4QKGL07sl5qXVS6cqW62w1QZbTXB6lyvESswgB6d3jXhp9dIJsNUKW22w1QSndymExGz6OL1rxEurl06ArVbYaoOtJji9SyEkZpDD1SSUowkzN61eOgG2WmGrDbaacHXMosRIzKbP1SSUowkzN61eOgG2WmGrDbaacHXMouQIfVJhoYvTuwF4afXSqcpWK2y1wVYTnN7lKvUKHlDo4vRuiXlp9dKpylYrbLXBVhOc3uUKsRIzyMHpXSNeWr10Amy1wlYbbDXB6V0KITGbPk7vGvHS6qUTYKsVttpgqwlO71IIiRnkcDUJ5WjCzE2rl06ArVbYaoOtJlwdsygxErPpczUJ5WjCzE2rl06ArVbYaoOtJlwdsyg5Qp9UWOji9G4AXlq9dKqy1QpbbbDVBKd3uUq9ggcUuji9W2JeWr10qrLVClttsNUEp3e5QqzEDHJweteIl1YvnQBbrbDVBltNcHqXQkjMpo/Tu0a8tHrpBNhqha022GqC07sUQmIGOVxNQjmaMHPT6qUTYKsVttpgqwlXxyxKjMRs+lxNQjmaMHPT6qUTYKsVttpgqwlXxyxKjtAnFRa6OL0bgJdWL52qbLXCVhtsNcHpXa5Sr+ABhS5O75aYl1YvnapstcJWG2w1weldrhArMYMcnN414qXVSyfAVitstcFWE5zepRASs+nj9K4RL61eOgG2WmGrDbaa4PQuhZCYQQ5Xk1COJszctHrpBNhqha022GrC1TGLEiMxmz5Xk1COJszctHrpBNhqha022GrC1TGLkiP0SYWFLk7vBuCl1UunKlutsNUGW01weper1Ct4QKGL07sl5qXVS6cqW62w1QZbTXB6lyvESswgB6d3jXhp9dIJsNUKW22w1QSndymExGz6OL1rxEurl06ArVbYaoOtJji9SyEkZpDD1SSUowkzN61eOgG2WmGrDbaacHXMosRIzKbP1SSUowkzN61eOgG2WmGrDbaacHXMouSwPGEQwBQAKwCsAjA7w+d7Ang4/vyfANR2dJu5pnfrP3a24ouVihug+GKl1n/s7KzXDa2x/rdaib0KtGol9mpj/W9DJ2XlpdVLpypbrbDVBltt1J9+s6LfagX2Kfqt1vrTb854PXCQg6tIy+6GgUoALwEYCaAawGIA9e2uMxPA3Pj30wA83NHtZtv0nTHt64pra9pMQuHaGj1j2tczXj+kxjOWKdDaZsAMaNXGM5aFTjuAl1YvnapstcJWG2y1ccaURxRV29u2Vm3XM6Y8csB1uenjKtYyG+QQkckAblTV98YfXxM/svjVtOv8Ir7OsyJSBWATgAGaIyrbIIfMGg4csvbAL3hjGPT2DNNcAVVJC/ah6oDLK9GCFj3w8pC8tHrpBNhqha022GpDDn4FeGvIgZ/o+wr0zbaXc5CDisXynL7BANalffxKfFnG66hqC4BtAA5vf0MiMkNEFojIgs2bN2f+bv3WFXZ5QPuQ+UTdbJeH5KXVSyfAVitstcFWI28NKuxyoiJwMcihqner6iRVnTRgwIDMV9o2tLDLA6rEga/NlOvykLy0eukE2GqFrTbYaqTvhsIuJyoCy03fegDpO64h8WUZrxM/vdsPwJbOfLMzNl0J7Gk3tbWnJrq8m5lxxioA7Z/B1vjy7sVLq5dOgK1W2GqDrTbOeNezQNXbbS+seju6nMiK1cmCAKoAvAxgBPYPcjS0u86/o+0gxyMd3W6u6d0zpn1dcdUwxQ2iuGpYtxziSGk8Y1nbCbNueKJxipdWL52qbLXCVhtstXHGlEcUfddF07t912Uc4lBVDnJwFW2ZviOHiLwfwLcQTfLeq6q3iMhN8V/g+SLSC8CDACYA2Apgmqq+nOs2sw1yEBERJREHOahYTMeZVPVnAH7W7rLr036/C8C5lg1ERERE5GSQg4iIiIi6hps+IiIiojLATR8RERFRGeCmj4iIiKgMcNNHREREVAa46SMiIiIqA9z0EREREZUBbvqIiIiIygA3fURERERlwPRt2CyIyGYAazq4Wn8Ar5Ugp9zwfi0+3qc2eL/a4P1qo6P7dbiqDihVDCWXu01fPkRkAd+nsPh4vxYf71MbvF9t8H61wfuVSoVP7xIRERGVAW76iIiIiMpAUjd9d4cOSCjer8XH+9QG71cbvF9t8H6lkkjkOX1ERERE1FZSH+kjIiIiojRuN30iMkVEVojIKhGZneHzPUXk4fjzfxKR2tJX+pPH/fpZEWkWkSUi8qSIDA/R6U1H92va9c4RERURTvLlIZ/7VUTOi//OLhWR75e60aM8/j8wTER+KyKL4v8XvD9Epycicq+IvCoiL2b5vIjInPg+XyIix5e6kZLP5aZPRCoBfAfA+wDUA/ioiNS3u9onAbyuqqMA3A7ga6Wt9CfP+3URgEmqeiyARwF8vbSV/uR5v0JE+gL4DIA/lbbQp3zuVxEZDeAaAO9S1QYAV5U81Jk8/75+AcAjqjoBwDQAd5S20qXvAZiS4/PvAzA6XjMA3FmCJiozLjd9AE4AsEpVX1bVPQB+CODsdtc5G8D98e8fBXCGiEgJGz3q8H5V1d+q6o74w+cADClxo0f5/H0FgJsR/eNkVynjHMvnfr0MwHdU9XUAUNVXS9zoUT73qwI4OP59PwAbStjnkqr+DsDWHFc5G8ADGnkOwCEiMrA0dVQuvG76BgNYl/bxK/FlGa+jqi0AtgE4vCR1fuVzv6b7JID/My1Khg7v1/ipnKGq+tNShjmXz9/XOgB1IvIHEXlORHI90kKRfO7XGwF8XEReAfAzAJ8qTVqiFfr/X6KCVYUOIJ9E5OMAJgE4NXSLdyJSAeC/AFwUOCWJqhA9XXYaokelfyci41T1jaBV/n0UwPdU9ZsiMhnAgyJyjKq2hg4jouy8PtK3HsDQtI+HxJdlvI6IVCF6CmJLSer8yud+hYicCeA6AFNVdXeJ2jzr6H7tC+AYAE+JyGoAJwGYz2GODuXz9/UVAPNVda+q/h3ASkSbQMoun/v1kwAeAQBVfRZAL0TvH0udl9f/f4m6wuum7y8ARovICBGpRnQi8fx215kP4BPx7z8C4DfKFyXsSIf3q4hMAHAXog0fz4/KT877VVW3qWp/Va1V1VpE50pOVdUFYXLdyOf/A48hepQPItIf0dO9L5cy0qF87te1AM4AABEZi2jTt7mklckzH8CF8RTvSQC2qerG0FGULC6f3lXVFhG5EsAvAFQCuFdVl4rITQAWqOp8AN9F9JTDKkQnz04LV+xDnvfrbQD6APhRPBezVlWnBot2IM/7lQqU5/36CwDvEZFmAPsAXK2qfMQ/hzzv1/8AME9EZiEa6riI/6jOTUR+gOgfIP3jcyFvANADAFR1LqJzI98PYBWAHQAuDlNKScZ35CAiIiIqA16f3iUiIiKiAnDTR0RERFQGuOkjIiIiKgPc9BERERGVAW76iIiIiMoAN31EREREZYCbPiIiIqIywE0fEeUkIu8UkSUi0ktEDhKRpSJyTOguIiIqDF+cmYg6JCJfRvRWW70BvKKqXw2cREREBeKmj4g6FL8H618A7ALwL6q6L3ASEREViE/vElE+Dkf0nst9ET3iR0REzvCRPiLqkIjMB/BDACMADFTVKwMnERFRgapCBxBR9yYiFwLYq6rfF5FKAH8UkX9V1d+EbiMiovzxkT4iIiKiMsBz+oiIiIjKADd9RERERGWAmz4iIiKiMsBNHxEREVEZ4KaPiIiIqAxw00dERERUBrjpIyIiIioD3PQRERERlYH/H9Jbr/VW/iLsAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "is4TRYp8Debh" }, "source": [ "## Discrete Initial and Boundary Conditions\n", "\n", "The discrete initial conditions are \n", "\\begin{equation}w[i,0]=x[i]\\times x[i], \\ \\ 0 \\leq x[i] \\leq 1\\end{equation}\n", "and the discrete boundary conditions are \n", "\\begin{equation}w[0,j]=t[j], w[10,j]=2-math.exp(t[j]), \\end{equation}\n", "where $w[i,j]$ is the numerical approximation of $U(x[i],t[j])$.\n", "\n", "The Figure below plots the inital and boundary conditions for $t[0]=0.$" ] }, { "cell_type": "code", "metadata": { "id": "8KcrR8R8Debh", "outputId": "b8d9307a-e64d-4b42-8f84-5a4c45dd74b0", "colab": { "base_uri": "https://localhost:8080/", "height": 304 } }, "source": [ "w=np.zeros((N+1,time_steps+1))\n", "b=np.zeros(N-1)\n", "# Initial Condition\n", "for i in range (1,N):\n", " w[i,0]=x[i]*x[i]\n", " \n", "\n", "# Boundary Condition\n", "for j in range (0,time_steps+1):\n", " w[0,j]=time[j]\n", " w[N,j]=2-math.exp(-time[j])\n", "\n", "fig = plt.figure(figsize=(8,4))\n", "plt.plot(x,w[:,0],'o:',label='Initial Condition')\n", "plt.plot(x[[0,N]],w[[0,N],0],'go',label='Boundary Condition t[0]=0')\n", "plt.xlim([-0.1,1.1])\n", "plt.ylim([-0.1,1.1])\n", "plt.title('Intitial and Boundary Condition',fontsize=24)\n", "plt.xlabel('x')\n", "plt.ylabel('w')\n", "plt.legend(loc='best')\n", "plt.show()" ], "execution_count": 3, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAEfCAYAAABRfP+PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU1f3/8dcne1jDvkMQEGRfogJuICpqFdG6Fhes1tZ+cWktFWpdq5WK1r1WbF2qVqn+ENGqKArigsouS0EBEQn7jpCEkJzfH/dOmISZZBIymUzyfj4e88jMXT9z5mY+c84991xzziEiIiLxJyHWAYiIiEjFKImLiIjEKSVxERGROKUkLiIiEqeUxEVEROKUkriIiEicUhKvIczM+Y/MCq6/1l9/SCXGdJe/zecra5uxZmbP++/prljHUhOY2RC/PNfGOhaJTLjvCjMb7U+fVcHt6lioACXxShL05T6rkrc7xE+GIyu4fl9//dGVGZeUT9AXXMlHoZntNbPFZvagmbWJdaxSucwsy8weMbNFZrbNzPLNbLuZfWFmD5hZz1jHGG3+8X+XmfWNdSw1jZJ49TcEuBMoK4mv9B/5Jab39dcfXcb6q/3195c7QimvzUGP7UBdoDdwC7DczAbFMDapJGZW38xeBeYCNwF9gAxgD9AQOB4YCyzxl4t3u/G+Q9aFmDca73uotCS+319/daVHVoMlxToAqRzOuW5HuP6wyopFSuecaxn82sxSgJ8AfweaAy+aWRen4RTjlpnVBz7BS9wHgWeAZ4EFzrlCM0sAegAXATcClwCXxijcSuGcewN44wjW/wo4ou+x2khJXCTGnHMHgDfMLB14GegEHAMsj2lgciQm4SXwXOA859z7wTOdc4XAErxa+KPAi1UfotQEak6vAmY2yz//OdrM0v1zQyvNLMfMtpjZq2bWpcQ6mWbm8JqgAK4KcT41M2j5kNOA5/yXp4RYf0jQsmE7tpnZyWb2qJl9aWYbzOyAH/d7ZnZh5ZRS5ewzuOOZmSWa2c3++eb9ZrbDzN42s6wy9n28mb3lL/+jfy7zJr/2FE1fBz2vGya2VDP7rV8uu/1jaKWZ/dXMWoZZp8wOhhamw17JzkZmdoJfhtv8fS82szFmZqVsu6F/vv87M8s1sx/M7Bkzaxu+KMDMmprZr83sTTNb4fcd2Gdmy/332zrMepmBY9x/PdDMXjezjWZWYN756Tv8ZeaVEcPV/nI/RPr5m9kADtWq/1gygZfknNuO1xITalv9zOwlf/95frlPN7OflrL/ov9lM2vsl9V3/vrZftm3KuM9jDLvnP2P/v/BR2YWMsagdQ7r2BaYBpziT3quxHfQ2qBly+zYZmZDzWyKmW3yvxM2mdkbZnZqKesUfTeaWXv//a/3y+M7/9hsUNp7q9acc3pUwgN4HnDArBDzZvnzbgQW+M9z8c4BOf+xHegUtE47YBPwoz8/x38d/GgXtHxgO5lB0zbhnadywIEQ6w8OWnatv9yQErHXC9q2wzuft7vEtKfDlMld/vzny1mWR7LPwOdwL/Be0HvfG7RuDjAozPqX4jV/BpbdidfPwAGvB23/rnK+p9GBbZayzGX+MgVAkxDzmwUdP4FjaE/Q6x3AwIp8DuHeF16fDOcfH6P9sikEdpX4PB4Js91WwLclyj7wWWwBrglsP8S6Dwatl4/3PxL82WwBeodYLzNomUuCPr9d/rHwCNDWL2cH9CqlXD4JHE/l+Kwnceh/Ou0IvlOuC4oxcCwGv/8XgcQQ6631518e9Hyff7wE1v0OaBRmv08ELVfg77eQQ99hgW2W/K4YTYnvQL/8N/nl7vD+j4O/g+aGOtbCxHVvUFyFJeJywP1h1gvMP8//TALfKflB8+YCyRX9rGL5iHkANeVBZEl8p//PMxxIxGsJOQn4wZ//nxDr3kUEiTDoYMwsMf2wf6ww64f7x6wDvIbXsa5x0PQM4P849IV8UUVjD7Hekewz8Dns9P9hLwZS/Hm98ZowHfBViHU74SUZB0wHjgqK57d4X6CB5HVXOd9T4HNwIeYlAyOAjaWVF/Auh5L1Rfhf4EAWXi3e4X0xNi3v50DZSXwfkAc8DrQI+jwe49CXao8Q233fn7/Vf48J/vSTgFVB5bk2xLo3AuOBXkCSPy0RGMChH2hLASuxXiaH/h/24v34yvTnJQU9f8df5q9hyqRL0HvrFK7sQqwX+NHyanmOkRLbGMyhBP4a0NafXg+4jUPJ64+l/C/vBBbi/2D13/sIf7oDHgix7qigspsIZPjTWwAv4CXjfUSYxIPmzfLnjS7lPQeOtVDHwqVBcT2Of4wDTYKOQQdcHmLd4B9BHwI9/empwM859OPm1xX9vGL5iHkANeVBZEl8P9A5xPyfcqhmlVJi3l3EMIlH8L6v8NebGWJeRLFXoKxL2+fzQWVxYoj5A4Lmty8x75/+9BWEqEEBfwxa965yxjw6aN3gmshWDn0h/wD8iRA1ArykF1h/eIj5LfCSuwPuKe/nQNlJ3AHPhFk38APijlJiHhpivc5BX6Bry1meqcAyf91TSszLDNrvp/g/HEJs43wO1ehDlfmfI/nfKbFOctC+xx/BMf5hUPyhatuB2PYCDUrMWxt0nIVq0bnFn7+mxHTD+2EV8ljx538Q9P6GhDnGDysvjiCJ+/sN/DB6Jcy6/+ZQC0NCiXmBeJcCqSHWfdyf/1FFP69YPnROvGq97pxbFWL6NLyDKBXviy2evOX/HWhmidVon5845z4tOdE5Nx9Y778suj7XP6d7gf/yYedcbohtPkLlXILXIujRFO9LCqAB0BhID7FOoB/APOfc9JIznXOb8Xq3g9f6EA33h5n+pv+35PXOgZi/cM7NLLmS/78wuSKBOOfy8BIKwAmlLPqQ8zqRhfIW3mV+zYBzgmf457+v9F8+W47QGgc931GO9YL33RgY6r+83zlXEGKxv+D9AKoHnB1mU5Ocd769pKn+345mFtz3oi9eaxSE+Kydl/H+XEb40dCXQ9+L94ZZ5m7/byZwXJhl/uofNyUFyiMur9dXEq9ac0NNdM7l49UGABpVXTiRMbMkM7vGvE5lG/0OIYGOQzv9xdKoxNgrYZ8hy9qX7f8NXvcovOZhgI9DreSc+xGYH9k7CM85Z8EPoDVeM+c64NfAx+ZdohSsv//3sGQY5CP/79Elvpwrww7n3Jow80KVJxyKOWR5RjAPM+tmZk+Y2ddmtse8wXECx8FN/mIhO7j55oSb4Zw7iNdEDHB1idnDgTZ4505fLy3GKOiH98POEf5Y3M2hY7F/qGUI/z+QHfQ8I+h5YDubnXMrw6z7Od5ppaoUiGurc25ZqAX8eLNLLF9SWeVR7b57I6FLzKrW3lLmBWp+yVURSKTMrB7e+eHBQZNzONQMDF6NErwe1duqyT7LW9bNgp5vKGXd7FLmVYhzbiPwlpnNBf6HV/O4ieK1jkB8pe0/0MJgeDX8fZUYZkWO3UDMFSpPM7sU+FfQdgvxOkYFalP18D7/0n6wbC1lHsA/gN8DZ5lZC79FA7xzpeCd1y5P60tw7btx2KVKFyi33f4Px3ACn3ezMPNDfmbOudygiwlC/Q+E/bycc3lmtg0IeSVElERy7INXHm0oZ3lw6PiNy3yomriU5Xa8ZLoNuAqvU1Md51xz5w1aEjxMaNjLjOJgnzHnnNuE12ELwjeJp1VRODFlZs3wBkhJxmtyz8Lrp9DIOdfSPw4eDiwebjthmqKD53+LV9tNwutrgZk1wWsZgfI1pQda1QKnzPqUZ90QUo9w/ZqmVhz75aUkLmW5yP97g3PuX865LSXmtyi5QpzuM7jGVlrzbGnzKkNgyMqjSkwPxNe+lHUD1107irdOBJo/S/sSbBhRdOUTiLki5XkWXk17OfAz59x8P0EGq6zj4B/+30CT+iggBVjmnPuyAtsLnPI43cwqkngC5Zbu/5gJJ/B5l9XaUN79hv28zBtdsGkl7S9SgbjalbFcZZdHXFASr/4CzccVrXEe6fqBf4yFYeafVsHtVrd9rsG73Ang5FAL+OeZSx0ophIEWhlKJqwF/t9TzMIOrBIY8OIb51xwU3rgfYUcXMXf3oDyBhqBQMwhy9N3SpjpgVi/DtUxzY857AAf5fQ6Xhl1N7PjOZTMn6vg9ib5fxvj9XEoU4nPdCHeDzE41MGt5PINOfSZLQi1TAUEttPCzI4Os8xgKtbsfCTfQ4G46ppZyE5rfrxtSixfKyiJV397/L8ZpS4VvfV3+397lZzhn7u+rYLbrVb79Hve/j//5c1mFqop80a8a8ajwswygNP9lyW/iAKdq3rgDVpRct0WwK/8l/8pMXuJ//fYMCN1jaLsWk5FvOb/HWRmhyVyMzsKbzCQUALHQM8wP1p+waGe1EfEvxLhJf/lQ3h9EvKp4FCozrl5HPoM7jOz00tb3m++/2/Q+js4VJu/NcxIcbfitaz8iHe9e2VYxKFTAbeGiNOAcRXc9pF8DwXH9Ycwy9zl/10LfFWBfcQtJfHqL9Ab80QrMTRrOdcP1DLKK3AZz1/NrKgWaGbH4l3L2qQC26yO+wTvsppcvHHLp5pZR3+/6WZ2M9413LtLWb9C/J74WXiXPDX3Jz8evIxz7hMOnS9/1swuDFxeZ94wn+/j9a7dDDxaYhef4XVWSgFeCXpfdczsl3jnnndSyfxL/AKf5etmdk4gIZnZCf77CXXJD8AMvNpoT+Ax/wcOZtbAzMYCT+IN5lNZAk3qgcvV3g5xGqc8foH34ykNeMfMnjSz/kHv38ysp5ndjXfXrrNKrH87Xu21P/Cq+UPUmlk9M/sDh5LpBOfcHiqB/0P2Lv/lz83sL0Hl3gKvf8CpVOwyy8D30AV+K0J54/qj//I8M3vc/+GDmTUxs8fwRjsEb/CbcJcU1kyxvlC9pjyIbLCX0aWsv5bQAygkc2gAhkK8S9HW+o+2QcuFHOzFn/dx0PztQesPjGD/R+GdYwqsn8OhoWD3A2eE2zcVH7HtSPYZ+BzuKmX7YT8Pyh529YWyth9mn6ODtlly+Nvg4R8LgbvDbKMZh5paA+VSctjVcMPJnk/xITx3B+33n+HKjTKGwizx3kId+yWHXd1P5MOu/jVovcBnEXgP73FoGM7nS6yXGVinnJ/R3KB9nVMJ3wkN8Fojgt/DQbz/wZKf+WH/I8Avg95vof/5Bh+bL1H6sKtDSomttO+LJ0rEu4MKDrsaNK8b3g8257/3bH87n0Z6rFF82NUCP67gY7qsYVcPe69HcrxUl4dq4tWc8zrzDMNr2svGq2118B+Rnpu6APgb3mhG9YLWL7PTjfOuDT4O7wtjC96wl7vw7rZ1rCvj5g4VEYt9Bu37Vbza2H/9fabgda66Ga/HuAu/dsRalHgcxDsn/yJwgnPuzjCxbQUGAb8D5uF9GabgJclH8IY9DXldtPNuE3kGXjPtXrwyXQRc45y7phLeU0jOu3zuWLyE/L2/3914Pxz6U8q9o51zv8UbP3whXgJI9J/fjHfDkMq+XnmK/3cj3hC3R8Q5t8c5dxHesfw4Xs18D15y3wN8AUwAjnHOjQ6x/tN4ZfdvP6Z6eGX3Ad6Qw5e7MnrfVzDuMXjjrn+JV+6GVxE4xzn3WAW3uQLvVNF7eO+hJd53UKk3wSmxjT/ifRe+iddxsx7eD6JpwGnOufEViS3emf9LRESkVjOzD/A6Tf7FOVfRc78iVUpJXERqPTPrDHzjvzzahR4eWaTaUXO6iNRq/hUPj+M1G7+tBC7xRDVxEamV/CsObsY7P5uKd2XCAOfc8pgGJlIOqomLSG2Vgde5qgDvxh5nKIFLvIm7mnjTpk1dZmZmrMMQERGpEvPnz9/mnAs5BG/c3bUlMzOTefPmxToMERGRKmFm34ebp+Z0ERGROKUkLiIiEqeUxEVEROJU3J0TDyU/P5/169eTm5sb61BEwkpLS6Nt27YkJyfHOhQRqSFqRBJfv3499evXJzMzk/C3WhaJHecc27dvZ/369XTs2DHW4YhIDVEjmtNzc3Np0qSJErhUW2ZGkyZN1FokIpWqRiRxQAlcqj0doyJS2WpMEhcREaltlMQrSb169cpc5tprr2X5cm9Uxz//+c/F5g0ePLjC+9i0aROXXnopnTp1YsCAAZx99tl88803IZctj1mzZnHOOecAMG3aNCZMmADA1KlTi94HwB133MGMGTOOeH8iIlI+NaJjW3lNXZjNxOkr2bArh9YZ6Ywd3pWR/dpEfb//+Mc/ip7/+c9/5g9/+EPR688//7xC23TOcf7553PVVVfx6quvArB48WI2b97M0UcffWQBBxkxYgQjRowAvCR+zjnn0L17dwDuueeeStuPiIhErtbVxKcuzGb8lCVk78rBAdm7chg/ZQlTF2ZXyvZnzZrFkCFDuPDCC+nWrRujRo0iMD79kCFDmDdvHuPGjSMnJ4e+ffsyatQo4FAt+8cff2TYsGH079+fXr168eabb5a6v5kzZ5KcnMyvfvWroml9+vThpJNOwjnH2LFj6dmzJ7169WLy5Mllxvjee+/RrVs3+vfvz5QpU4q2+fzzzzNmzBg+//xzpk2bxtixY+nbty+rV69m9OjRvP766wB8+OGH9OvXj169evHzn/+cvLw8wBsu98477yx6XytWrKiM4hYRqdVqZBK/5Ok5vDbvBwDyCwq55Ok5vLFwPQAPvLeCnPyCYsvn5Bdwz9vLANix7wCXPD2HGcs3A7Blb/l7Ey9cuJBHHnmE5cuXs2bNGj777LNi8ydMmEB6ejqLFi3i5ZdfLjYvLS2NN954gwULFjBz5kxuueUWSrtJzdKlSxkwYEDIeVOmTGHRokUsXryYGTNmMHbsWDZu3Bg2xtzcXH7xi1/w1ltvMX/+fDZt2nTYNgcPHsyIESOYOHEiixYtolOnTkXzcnNzGT16NJMnT2bJkiUcPHiQp556qmh+06ZNWbBgAddffz0PPvhg2QUpIiKlqpFJvDQbd4dOyjv25VfaPo477jjatm1LQkICffv2Ze3atRGv65zjD3/4A7179+a0004jOzubzZs3VyiOTz/9lMsuu4zExERatGjBKaecwty5c8PGuGLFCjp27EiXLl0wMy6//PJy7W/lypV07NixqBn/qquuYvbs2UXzL7jgAgAGDBhQrjIREZHQonZO3MyeBc4BtjjneoaYb8CjwNnAfmC0c25BZex78i8HFT1PTkwo9rp1RjrZu3IOW6dNRjoAjeumFFu+ef20cu8/NTW16HliYiIHDx6MeN2XX36ZrVu3Mn/+fJKTk8nMzCz12uIePXoUNWVXVYwVFdhnVe1PRCRSv3nrCZ5a8Cfy3FZSrRnX97+dh88dE+uwyhTNmvjzwJmlzD8L6OI/rgOeKmXZSjN2eFfSkxOLTUtPTmTs8K5VsfsiycnJ5OcfXvvfvXs3zZs3Jzk5mZkzZ/L992HvQAfAqaeeSl5eHpMmTSqa9vXXX/PJJ59w0kknMXnyZAoKCti6dSuzZ8/muOOOC7utbt26sXbtWlavXg3AK6+8EnK5+vXrs3fv3sOmd+3albVr17Jq1SoAXnzxRU455ZRS4xcRibXfvPUEj87/HXlsAXPksYVH5/+O37z1RKxDK1PUkrhzbjawo5RFzgP+5TxfABlm1ipa8QSM7NeG+y/oRZuMdAyvBn7/Bb2qpHd6sOuuu47evXsXdWwLGDVqFPPmzaNXr17861//olu3bqVux8x44403mDFjBp06daJHjx6MHz+eli1bcv7559O7d2/69OnDqaeeygMPPEDLli3DbistLY1Jkybxk5/8hP79+9O8efOQy1166aVMnDiRfv36FSX8wPrPPfccF110Eb169SIhIaFYhzsRkeroqQV/wllesWnO8nhqwZ9iFFHkrLROU0e8cbNM4O0wzelvAxOcc5/6rz8EbnXOzQux7HV4tXXat28/oGTt9H//+x/HHHNMpccvUtl0rIpUP3ZXAliIXOgMd1dh1QdUgpnNd85lhZoXFx3bnHOTnHNZzrmsZs2axTocERGpQVIsdF5JDTO9OollEs8G2gW9butPExERqTI/7zUec6nFpplL5fr+t8coosjFMolPA640z0Bgt3NuYwzjERGRWuSHHfspLHQ8dcHN3NR/Iqk0B2ek0pybBjwYF73To3mJ2SvAEKCpma0H7gSSAZxzfwfewbu8bBXeJWZXRysWERGRYOt37ufsRz/hupOP4oZhXXh4xA08POKGWIdVblFL4s65y8qY74D/i9b+RUREwmmTkc71Qzsxsm/VXplU2eKiY5uIiEhleGfJRjbsysHM+PWQzrT2B/qKV0rilSQxMZG+ffvSp08f+vfvX+G7kpVHZmYm27Ztq/Tt5ufnM27cOLp06UL//v0ZNGgQ7777bqVsO3Cjlw0bNnDhhRcCsGjRIt55552iZYJve3okdu3axd/+9rei17NmzaJhw4acffbZRdNeeOEFunTpQpcuXXjhhReKpg8dOpR69eoxb95hVzyGtGPHDk4//XS6dOnC6aefzs6dO484fhGpXLv2H+DW//c1j3/0baxDqTS1Mom/vORlMh/JJOHuBDIfyeTlJS+XvVIZAjc0Wbx4Mffffz/jx4+vhEgrV0FBQdkLAbfffjsbN25k6dKlLFiwgKlTp4Ycoe1ItG7dumi42JJJfMSIEYwbN+6I91EyiQOcdNJJRfvasWMHd999N19++SVfffUVd999d1HynTlzJllZIS/LDGnChAkMGzaMb7/9lmHDhlXKjxARqVwZdVKYfN0g7jy3R6xDqTS1Lom/vORlrnvrOr7f/T0Ox/e7v+e6t66rlEQesGfPHho1agRQ6u1AzznnnKJ1xowZw/PPPw+Ev23n9u3bOeOMM+jRowfXXnttsbubjRw5kgEDBtCjR49iQ7DWq1ePW265hT59+nDfffcxcuTIonkffPAB559/frHY9+/fzzPPPMPjjz9eNNZ5ixYtuPjiiwFvKNZevXrRs2dPbr311mL7ue222+jTpw8DBw4sumnLd999x6BBg+jVqxd//OMfi5Zfu3YtPXv25MCBA9xxxx1MnjyZvn37Mnny5KLbngaWO/XUU+nduzfDhg1j3bp1AIwePZobb7yRwYMHc9RRR4UcP37cuHGsXr2avn37Mnbs2MPmT58+ndNPP53GjRvTqFEjTj/9dN57770Qn2jZ3nzzTa666irAu/HL1KlTK7QdEalchYWOu6YtK7rddPfWDUgrMfR2PKt1Sfy2D29jf/7+YtP25+/ntg9vO6LtBu4P3q1bN6699lpuv927vrC024GWJtRtO++++25OPPFEli1bxvnnn1+U0ACeffZZ5s+fz7x583jsscfYvn07APv27eP4449n8eLF3H777axYsYKtW7cC8Nxzz/Hzn/+82H5XrVpF+/btadCgwWExbdiwgVtvvZWPPvqIRYsWMXfu3KJktW/fPgYOHMjixYs5+eSTeeaZZwC46aabuP7661myZAmtWh0+qm5KSgr33HMPl1xyCYsWLeKSSy4pNv+GG27gqquu4uuvv2bUqFHceOONRfM2btzIp59+yttvvx2y5j5hwgQ6derEokWLmDhx4mHzs7Ozadfu0FAFbdu2JTs79FAFJ510En379j3sMWPGDAA2b95c9P5atmxZ4TvPiUjlyi8sZOWmvfxv055YhxIVtS6Jr9u9rlzTIxVoTl+xYgXvvfceV155Jc65Um8HWppQt+2cPXt20e1Bf/KTnxTV9gEee+yxolrwDz/8wLffeud8EhMT+elPfwp446xfccUVvPTSS+zatYs5c+Zw1llnRfwe586dy5AhQ2jWrBlJSUmMGjWq6FajKSkpRS0LwTF/9tlnXHaZd6HCFVdcEfG+AubMmcPPfvazovU//fTTonkjR44kISGB7t27Rz1pfvLJJyxatOiwx2mnnXbYsmaGd5M+EYmVnAMF5OYXkJqUyPM/P5ZxZ5Z+H4p4FbVLzKqr9g3b8/3uw+8M1r5h+0rbx6BBg9i2bVtRjTeUpKQkCgsPjclb8naj5blt56xZs5gxYwZz5syhTp06DBkypGh7aWlpJCYeajq6+uqrOffcc0lLS+Oiiy4iKan4IdC5c2fWrVvHnj17QtbGw0lOTi5KXCVjjlZCC76dakXuAdCmTRtmzZpV9Hr9+vUMGTIk5LInnXRSyH4BDz74IKeddhotWrRg48aNtGrVio0bN4a9eYyIRF9BoeOq574iIz2Zp68YQGpSzWk+L6nW1cTvG3YfdZLrFJtWJ7kO9w27r9L2sWLFCgoKCmjSpEnY24F26NCB5cuXk5eXx65du/jwww/L3O7JJ5/Mv//9bwDefffdok5Yu3fvplGjRtSpU4cVK1bwxRdfhN1G69atad26Nffeey9XX334+Dp16tThmmuu4aabbuLAgQMAbN26lddee43jjjuOjz/+mG3btlFQUMArr7xS5q1GTzjhBF599VXAu1d6KOFubQowePDgYuufdNJJpe4v0u0CDB8+nPfff5+dO3eyc+dO3n//fYYPHx5y2bJq4iNGjCjq3f7CCy9w3nnnRRyniFSuxATj3D6tGdG3dY1vFat1SXxUr1FMOncSHRp2wDA6NOzApHMnMarXqLJXLkXgnHjfvn255JJLeOGFF0hMTAx7O9B27dpx8cUX07NnTy6++GL69etX5j7uvPNOZs+eTY8ePZgyZQrt23utB2eeeSYHDx7kmGOOYdy4cQwcOLD0Mhg1inbt2oW9m9a9995Ls2bN6N69Oz179uScc86hQYMGtGrVigkTJjB06FD69OnDgAEDykxWjz76KE8++SS9evUKe7556NChLF++vKhjW7DHH3+c5557jt69e/Piiy/y6KOPlrq/YE2aNOGEE06gZ8+eITu2NW7cmNtvv51jjz2WY489ljvuuIPGjRtHvP1g48aN44MPPqBLly7MmDGjUnrXi0j5/LBjP0uzdwNwxcAOnNO7dYwjir6o3oo0GrKyslzJa3d1e8fyGTNmDP369eOaa66JdShVatasWTz44IO8/fbbES0/ZMgQHnzwwXJdalYWHasi0eGc48K/z2Hn/gN88JtTSEyoOTXw0m5FWuvOidd2AwYMoG7dujz00KLZTdsAACAASURBVEOxDqXKpaSksHTpUs4+++xi16WHMnToUNasWUNycnIVRSciR8LMePCiPhw4WFijEnhZlMRrmfnz58c6hJgZPHhwUa/5ssycOTO6wYhIpfjk2618vX43/ze0Mx2b1o11OFWuxpwTj7fTAlL76BgVqXzvLt3EW4s3kJsf2YiUNU2NqImnpaWxfft2mjRpUuN7Ikp8cs6xfft20tLSYh2KSI2Qd9C7BvyeET3IyS+oUaOwlUeNSOJt27Zl/fr1pV6XLRJraWlptG3bNtZhiMS9f3yyhjcWZvPqdQOpn5ZM/cQa06hcbjUiiScnJ9OxY8dYhyEiIlWgc/N6dGtZs8ZAr6gakcRFRKRmO1hQyPKNe+jdNoMhXZszpKtGRYQa1LFNRERqrkc//JYLn5rDDzv2l71wLaKauIiIVHvXnNiRDk3q0q5xnbIXrkVUExcRkWppx74DPPzBNxQUOjLqpHDhAHUMLUlJXEREqqXpyzbx949Xs3JT+BsZ1XZqThcRkWrFOYeZcdlx7TmhU1PaN1ETejiqiYuISLWxNHs3P3ns06IObErgpVMSFxGRaiM5MQEzOFioYYojoeZ0ERGJuXXb99O+SR26tqzP2zecqCG0I6SauIiIxNSc1ds59aFZvLd0E4ASeDmoJi4iIlVq6sJsJk5fyYZdObTOSOc3p3Xh10M6cWKXprEOLe6oJi4iIlVm6sJsxk9ZQvauHByQvSuH299cxlHN6lEvVfXK8opqEjezM81spZmtMrNxIea3N7OZZrbQzL42s7OjGY+IiMTWxOkrySlx7++c/AImTl8Zo4jiW9SSuJklAk8CZwHdgcvMrHuJxf4I/Mc51w+4FPhbtOIREZHY27Arp1zTpXTRrIkfB6xyzq1xzh0AXgXOK7GMAxr4zxsCG6IYj4iIxNB32/aRmBC601rrjPQqjqZmiGYSbwP8EPR6vT8t2F3A5Wa2HngHuCHUhszsOjObZ2bztm7dGo1YRUQkyhqmJ9O2UTqpScVTT3pyImOHd41RVPEt1h3bLgOed861Bc4GXjSzw2Jyzk1yzmU557KaNWtW5UGKiEjFOOeYvmwThYWOxnVTmPm7Ifzlp71pk5GOAW0y0rn/gl6M7FeyjieRiGZXwGygXdDrtv60YNcAZwI45+aYWRrQFNgSxbhERKSKzPpmK798cT6PX9aPc/u0xswY2a+NknYliWZNfC7Qxcw6mlkKXse1aSWWWQcMAzCzY4A0QO3lIiJxzjlv2NQhRzfj6SsG8JNerWIcUc0UtSTunDsIjAGmA//D64W+zMzuMbMR/mK3AL8ws8XAK8BoF/jkRUQkLi3fsIcLnvqcTbtzMTOG92hJQpgObXJkonplvXPuHbwOa8HT7gh6vhw4IZoxiIhI1UpIgL25B9m5/wAtG6bFOpwaLdYd20REpAZwzvHVdzsA6NayAdNvPpljWjUoYy05UkriIiJyxF7+ch0XPz2Hhet2AoS9HlwqlwaqFRGRI3bhgLakJiXQt11GrEOpVVQTFxGRCln8wy6u+9c8cvMLSEtO5KKsdrqNaBVTEhcRkQrZtCeXFZv2smVPXqxDqbXUnC4iIhErLHR8u+VHurasz/AeLRnStRmpSYmxDqvWUk1cREQi9vCMbxj55Gdk+3cdUwKPLdXERUQkYlcM7EDLhmm01vXf1YJq4iIiUqp5a3dwz1vLcc7RvEEao47voA5s1YSSuIiIlGrO6u18tGIzO/fnxzoUKUHN6SIicpiCQsfmPbm0zkhnzKmdGX1CJvXTkmMdlpSgmriIiBxm/JSvuejvc/gx7yBmpgReTakmLiIih7l8YAf6tMugXqrSRHWmT0dERAD4fPU2Vm/5kSsGZdK7bQa922oI1epOzekiIgLAq1/9wEtfrCPvYEGsQ5EIqSYuIlKL5RcUsv9AAQ3Tk5nw014UFDoN4BJHlMRFRGop5xy/enE+e/MO8sovBlInRSkh3ugTExGppcyM8/q1IS+/QPf/jlNK4iIitcyslVsAGNK1OSP6tI5xNHIklMRFRGqRwkLHQ+9/Q3pKIqcc3UzDp8Y5JXERkRps6sJsJk5fyYZdObTKSOP3w7vxj6uyqJeapAReA+gSMxGRGmrqwmzGT1lC9q4cHLBhVy7jpyxhzurt1NUgLjWCkriISA01cfpKcvKLX/Odk1/AxOkrYxSRVDYlcRGRGqig0JG9KyfkvA1hpkv8URIXEamBdu4/QLirxlpnpFdtMBI1SuIiIjXI4h924Zyjab1U7jy3O+nJxb/m05MTGTu8a4yik8qmJC4iUkN8tmob5z35GW9/vRGAqwZ35P4LetMmIx0D2mSkc/8FvRjZr01sA5VKY8656G3c7EzgUSAR+IdzbkKIZS4G7gIcsNg597PStpmVleXmzZsXhWhFROJTQaEjMcEoLHS89OX3XHJsO41/XoOY2XznXFaoeVGriZtZIvAkcBbQHbjMzLqXWKYLMB44wTnXA7g5WvGIiNRE7y7ZyNmPfsLunHwSEowrB2Uqgdci0WxOPw5Y5Zxb45w7ALwKnFdimV8ATzrndgI457ZEMR4RkRqnRcM0mtVP1e1Da6loJvE2wA9Br9f704IdDRxtZp+Z2Rd+8/thzOw6M5tnZvO2bt0apXBFROLD1+t38epX6wDo374RL117PM3rp8U4KomFWHdsSwK6AEOAy4BnzCyj5ELOuUnOuSznXFazZs2qOEQRkerl2U+/42+zVpObr9p3bRfNcfeygXZBr9v604KtB750zuUD35nZN3hJfW4U4xIRiTvbfszDOWhWP5V7RvaksNCRlqxz37VdNGvic4EuZtbRzFKAS4FpJZaZilcLx8ya4jWvr4liTCIicSe/oJAL/vY54/7f1wA0SEsmo05KjKOS6iBqNXHn3EEzGwNMx7vE7Fnn3DIzuweY55yb5s87w8yWAwXAWOfc9mjFJCIST5xzmBnJiQmMO6sbRzWrG+uQpJqJ6nXi0aDrxEWkNtiyN5dfvTifG4d1YUjX5rEOR2IoJteJi4hIxTVISyYpIYG8g4WxDkWqMSVxEZFqIudAAU989C0HDhaSlpzI5F8OZHiPlrEOS6oxJXERkWriy++289AH3/DZqm0AmIW5DZmIT0lcRCSGnHN8t20fAEO6NmfGb09haDedA5fIKImLiMTQox9+yzmPfUL2rhwAOjWrF+OIJJ5Ec7AXEREJI3D52EVZ7WiQlkyrBho2VcpPNXERkSr25MxVjJ+yBPDu8f3zEzuSkKDz31J+SuIiIlUsL7+AnPwC8gt0+ZgcGTWni4hUgRnLN9OyYRo92zTk5tOOxky9z+XIqSYuIhJlufkF/HHqUp76eDUACQmmBC6VIqKauJm9BHwMfOKcWxHdkEREaobsXTm0bphGWnIiL117PO0ap8c6JKlhIm1O/ydwEvC4mXUCFgKznXOPRi0yEZE4M3VhNhOnr2TDrhya1U9l5/4D3H5Od64clEnn5rp0TCpfREncOTfTzGYDxwJDgV8BPQAlcRERvAQ+fsoScvILANiyN4+kBKMwzm4yJfElonPiZvYh8BlwCbASONY51y2agYmIxJOJ01cWJfCAg4WOZ2Z/F6OIpDaItGPb18ABoCfQG+hpZjq5IyLi2+CPuBbpdJHKEGlz+m8AzKw+MBp4DmgJpEYtMhGROPDmomx+2LGf1hnpRUOnBmudofqORE+kzeljzGwyXoe284BngbOiGZiISDz4Ys12Zq3cyi2ndyE9ObHYvPTkRMYO7xqjyKQ2iLR3ehrwV2C+c+5gFOMREanWDhYU8uxn33FG95ZkNq3LHef0ICUpgcQEIyEhoah3euuMdMYO78rIfm1iHbLUYJE2pz8Y7UBEROLBjn0HePzDVeQcKOSm07qQnnKo9j2yXxslbalSGnZVRKQM+/IO8s6SjVyU1Y7mDdJ49+aTaKNz3VINaNhVEZEyvPLVOsa+/jUrNu0BoG2jOho2VaoF1cRFRELYsjeXnfvy6dqyPlcOymRAh0Z0a9kg1mGJFKMkLiJSgnOOq5+bC8DbN5xISlIC/do3inFUIodTEhcR8X2/fR9tG9UhMcG457weZNRJUbO5VGs6Jy4iAqzctJfT/zqbf3/5PQADOjSmUzPdtESqNyVxEanVdufkA3B0i3rcfHoXzuzZKsYRiUROSVxEaq1/fLKGYQ99zK79BzAzfj2kM83qazRpiR9RTeJmdqaZrTSzVWY2rpTlfmpmzsyyohmPiIhzjoMFhQCc0LkpI/u2JiVJ9RmJT1E7cs0sEXgSb4z17sBlZtY9xHL1gZuAL6MVi4gIeEOmXvfifB6YvhKAY1o14I/ndKdOivr4SnyK5s/P44BVzrk1zrkDwKt4N08p6U/AX4DcKMYiIkJSYgLtGtWhuZrMpYaIZhJvA/wQ9Hq9P62ImfUH2jnn/lvahszsOjObZ2bztm7dWvmRikiN9d22fVw6aQ7fbdsHwB3ndufak46KcVQilSNmJ4LMLAHvzmi3lLWsc26Scy7LOZfVrFmz6AcnIjVGnZRENu/JI3vn4ff6Fol30Uzi2UC7oNdt/WkB9YGewCwzWwsMBKapc5uIHKnPV23jvv8uB6BFgzRm/PYUTuzSNMZRiVS+aCbxuUAXM+toZinApcC0wEzn3G7nXFPnXKZzLhP4AhjhnJsXxZhEpBZY+MMuZvxvC7v3e9eAJyZo1DWpmaLWJdM5d9DMxgDTgUTgWefcMjO7B5jnnJtW+hZEREKbujCbidNXsmFXDq0z0vndGUeDQbtGdcjKbMx1Jx/FNSd2JC05seyNicQxc87FOoZyycrKcvPmqbIuUltNXZjN+ClLyMkvKJqWlpxAWlIip3Zrzl8v6RvD6EQqn5nNd86FPNWsEQ5EJK5MnL6yWAIHyM0vJDU5gYkX9YlRVCKxoSQuInFlw67Qvcy37MnTuW+pdZTERSQuOOf4+JutNKqbHHJ+64z0Ko5IJPaUxEUkbjzx0bfUT00iPbn4V1d6ciJjh3eNUVQisaMkLiLV1vfb9/G71xazNzcfM+Pxy/rzwW+HcP8FvWmTkY4BbTLSuf+CXozs16bM7YnUNBr1X0Sqrd05+by3dBMX9G/D4E5NadkwDYCR/dooaYugJC4i1Yhzjofe/4akROPm046md9sM5ow/lfppoc+Di9R2SuIiEnOFhY6EBMPM2LA7h5TEQ2f6lMBFwtM5cRGJqa++28HQh2axfud+ACZe2IcJP+0d46hE4oOSuIhUOecc+/IOAtCucTotGqTxo/9a13qLRE7N6SJSpZxzXPHPr8iok8wTP+tPq4bp/OeXg2IdlkhcUhIXkSqxfud+2jaqg5kxtFtz6qbo5iQiR0rN6SISddOXbeLkB2ayYN1OAK45sSOXHtc+xlGJxD8lcRGJiu0/5rF6648AnNi5KWNO7cJRTevGOCqRmkXN6SJS6ZxzXDLpCxrVSea1Xw2mbmoSvz396FiHJVLjKImLSKU4cLCQd5ZsZESf1iQkGHed24OWDVNjHZZIjaYkLiKV4v3lm7h58iIa103h5KObcWKXprEOSaTGUxIXkQr75NutHCxwDO3WnLN6tuLf16YwqFOTWIclUmsoiYtIqaYuzGbi9JVs2JVD64x0xg7vysh+bXDO8cB7K6mbmsjQbs1JTDAGd1btW6QqKYmLSFhTF2YzfsoScvILAMjelcMtry3mwMFCLj62HU9d3p9m9XXeWyRWdImZiIQ1cfrKogQeUFDoeGD6CgDaNqpDapIGbRGJFSVxEQlrw66ckNO3/3igiiMRkVCUxEWkmPyCQhb6I6u1zkgPuUy46SJStZTERaSYidNXcsmkL9j+Yx5jh3clPbl4c3l6ciJjh3eNUXQiEkwd20Rquc17cnlkxrdce1JHOjWrx6jj2zPwqMY0qpPCyH5tAEL2TheR2FMSF6mFnHPsyT1Iw/RkEsx4++sNHN+xMZ2a1aNDk7p0aHJojPOR/dooaYtUU0riIrWMc45Lnv6CRnWTefqKLJrVT2XubaeRlqxe5iLxJqrnxM3sTDNbaWarzGxciPm/NbPlZva1mX1oZh2iGY9IbbVpdy7Pf/YdAGbGyH5tOLtXq6L5SuAi8SlqSdzMEoEngbOA7sBlZta9xGILgSznXG/gdeCBaMUjUhs55wBvXPN73l7OGv/WoD87vj3n9VUTuUi8i2ZN/DhglXNujXPuAPAqcF7wAs65mc65/f7LL4C2UYxHpNbYvCeXC/72GdOXbQbgwgFt+XjsUI5qVi/GkYlIZYpmEm8D/BD0er0/LZxrgHdDzTCz68xsnpnN27p1ayWGKFJz7N6fz9Ls3QA0rZdKesqhJvI6KUm0a1wnVqGJSJRUi45tZnY5kAWcEmq+c24SMAkgKyvLVWFoInHjVy/NZ/OeXD685RQSE4yXrx0Y65BEJMqimcSzgXZBr9v604oxs9OA24BTnHN5UYxHpEZZmr2bp2evYeKFvUlLTmTsmV1JS0rEzGIdmohUkWg2p88FuphZRzNLAS4FpgUvYGb9gKeBEc65LVGMRaRGyM0vYF/eQQB+zDvInNXbWLXF66zWv30jurduEMvwRKSKRS2JO+cOAmOA6cD/gP8455aZ2T1mNsJfbCJQD3jNzBaZ2bQwmxOp9Xbn5HPChI/456fepWLHd2zMZ+NOpWebhjGOTERiJarnxJ1z7wDvlJh2R9Dz06K5f5F4M3VhdrEhTi89rh2tGqZz4YC2NExP5qrBmZzQuQngXe+t24CK1G7VomObiHgJfPyUJUX3787elcMjM74lPTmBEX1ak5KUwI3DusQ4ShGpTnQXM5FqYuL0lUUJPKCg0FE/LZmUJP2risjhVBMXiaHd+/N5Yua3nHZMCzbsygm5zKbduVUclYjEC/28F6liyzfsYf73OwFIS0lgyoJslm3YQ+uM9JDLh5suIqKauEgV2LnvAI3qpgDwm8mLaFgnmf/8chCpSYl8Pv5UUpMSaVw3pdg5cYD05ETGDu8aq7BFpJpTEheJsrvfWsY7SzYyZ9wwEhKMhy7uQ6uGaUXzAz3MA/fsDu6dPnZ4V93LW0TCUhIXqWRzVm/nL++t4PmrjyWjTgqnHdOCDo3rkF9YSGpCYqnXdY/s10ZJW0QipnPiIkdoT24+z376Hau27AWgbmoizjk27/FGET6hc1NGn9BR13SLSKVTTVykArb9mMe+vIN0aFKXgwWO+975H4WuG52b16d32wzeHHNirEMUkVpASVwkQgWFjsQEwznHTx77hAEdGvG3UQNoXDeFT28dSquG6kUuIlVLSVwkAve/+z++WLODN//vBMyM+0b2KnZ/biVwEYkFJXGptUqOUx7cE3zBup386/O1PHBhH1KSEji6eX2cg4MFhSQlJnBa9xYxjl5ERElcaqlQ45T//vXF7MnN58pBmezcd4DPV29n3Y59dG5en58OaBvjiEVEDqfe6VIrhRqn/ECB46H3vwFgSNfmzBk/jM7N68ciPBGRiCiJS62xZU9u0WVg4cYp35OTD0BigpGYYFUWm4hIRSiJS42VvSuHeWt3FL2+4p9fcfdby4Hw45FrnHIRiSdK4lJjrN+5n+nLNhW9vvft5dz06qKi17ef051bz+wGwNjhXUlPLj74isYpF5F4o45tEreyd+Uwa+UWLju2PQkJxitfrePvH6/h6zvPoG5qEjec2gUA5xxmxoldmhatq3HKRaQmUBKXuLF5Ty5vLd7A+f3a0KReKl+s3s5tbywlq0Njurasz8+O78D5/dpSJ8WrYXdv3aDU7WmcchGJd2pOl2pr1/4D/G3WKlZs2gN4ndHu/e//WLBuFwCn92jB7LFDObpFPQDaZKTTuXk9zNQhTURqB9XEpUqUNrBKQN7BAv4+aw292zVkaNfmADw4fSX1UpPo1rIBvdo0ZM74U4tGR2uQlkyDtOQqfy8iItWFkrhEXaiBVcZPWQLAuh37aVw3hcsHdiAlMYGXv/yevINtGdq1ORl1Ulhw++lk1EkBICkxQcObiogEURKXqHtg+orDBlbJyS9g4vSVdGhShzYZ6Vw+sANmxuzfDyUtqNd4IIGLiMjhlMSlUi3bsJvlG/ZwUVY7AO5+axkbduWGXHbDrhw++f1QEoIGVUlL1j23RUQipY5tUi479x3g81XbKCx0ALw+fz3DH57NwYJCAN5dsolxU5aQ778+LrMx9dNC/1ZsnZFeLIGLiEj5KIlLMc45tu7NI++g1/y9YN1Orn9pPlv35gHwztKN/OwfX7Jxj1e7rp+WRMemddmX5y1/9QmZfDF+GEl+cj6rVyv+dF5PDawiIhIFSuI1zNSF2Zww4SM6jvsvJ0z4iKkLs0Mu55xXk964O4e/fvAN323bB8Csb7Zy7H0zWJq9G4DcAwWs3LyXHfsOAHBqt+a8fO3xNKnrnase3qMlf79iAA3reL3Em9RLpVn91GKXeY3s14b7L+hFm4x0DO9SsPsv6KVrtEVEjpAFvszjRVZWlps3b16sw6iWSvYCB0hLTuCSrHZcMagDnZvX5/vt+xj55Gfcc15Pzu3TmlVb9nL6w7N5alR/zuzZii17c/nv1xs5s2dL9QQXEakGzGy+cy4r1Lyo1sTN7EwzW2lmq8xsXIj5qWY22Z//pZllRjOeyvabt54g7e4W2F0JpN3dgt+89USlbHdvbj57cvOLXn/13Q6+3by36PVzn33HZ6u2Fb3+7eRFvLkoO+TtNXPzC3lhzve8v3wzAM3rp3Fmz1ZFN/o4qmk9VvzpTM7s2apo/tUndFQCFxGJA1FL4maWCDwJnAV0By4zs+4lFrsG2Omc6ww8DPwlWvFUtt+89QSPzv8deWwBc+SxhUfn/47fvPUE87/fUezuWZPnrivWrH33W8t4cuaqoteXTprDHW8uLXp9xsOz+ZN/ty2AX788n2c/W1v0+tEPv+X9oBt9LN+4h027c8PeXtOAqwd3BCA9JZH7L+jFgA6NAEhIMFKT1CNcRCQeRfMSs+OAVc65NQBm9ipwHrA8aJnzgLv8568DT5iZuTho439qwZ9wlldsmrM8nlrwJ7KzB5CQAK9eNwiAV+f+QP205KJzwFv25GEcOmfcr30j2gTdAvO3px9drCY86cosmtZNLXr96a2nFuso9t7NJwPwrznfkx0ikbfOSCc9RYlaRKSmido5cTO7EDjTOXet//oK4Hjn3JigZZb6y6z3X6/2l9lWYlvXAdcBtG/ffsD3338flZjLw+5KAAtRds749te7SUxIoGPTugDkFxSSlGBRH9M71Dnx9OREdSITEYljMTsnXlmcc5Occ1nOuaxmzZrFOhwAUi10HKnWjM7N6xclcIDkxIQquSmHeoGLiNQu0WxOzwbaBb1u608Ltcx6M0sCGgLboxhTpbm+/+08Ov93xZrUzaVy/YDbYxiVbq8pIlKbRLMmPhfoYmYdzSwFuBSYVmKZacBV/vMLgY/i4Xw4wMPnjuGmAQ+SSnNwRirNuWnAgzx87piyVxYREakEUb1O3MzOBh4BEoFnnXP3mdk9wDzn3DQzSwNeBPoBO4BLAx3hwtF14iIiUpuUdk48qjdAcc69A7xTYtodQc9zgYuiGYOIiEhNFRcd20RERORwSuIiIiJxSklcREQkTimJi4iIxCklcRERkTilJC4iIhKnlMRFRETilJK4iIhInFISFxERiVNK4iIiInFKSVxERCROKYmLiIjEqajexSwazGwr8H2s4yihKbAt1kHECZVVZFROkVE5RU5lFZnqWE4dnHPNQs2IuyReHZnZvHC3iZPiVFaRUTlFRuUUOZVVZOKtnNScLiIiEqeUxEVEROKUknjlmBTrAOKIyioyKqfIqJwip7KKTFyVk86Ji4iIxCnVxEVEROKUkriIiEicUhIvBzM708xWmtkqMxsXYn6qmU32539pZplVH2XsRVBOvzWz5Wb2tZl9aGYdYhFndVBWWQUt91Mzc2YWN5e+VKZIysnMLvaPq2Vm9u+qjrG6iOD/r72ZzTSzhf7/4NmxiDOWzOxZM9tiZkvDzDcze8wvw6/NrH9Vxxgx55weETyARGA1cBSQAiwGupdY5tfA3/3nlwKTYx13NS2noUAd//n1tbGcIi0rf7n6wGzgCyAr1nFXx3ICugALgUb+6+axjrsal9Uk4Hr/eXdgbazjjkE5nQz0B5aGmX828C5gwEDgy1jHHO6hmnjkjgNWOefWOOcOAK8C55VY5jzgBf/568AwM7MqjLE6KLOcnHMznXP7/ZdfAG2rOMbqIpJjCuBPwF+A3KoMrhqJpJx+ATzpnNsJ4JzbUsUxVheRlJUDGvjPGwIbqjC+asE5NxvYUcoi5wH/cp4vgAwza1U10ZWPknjk2gA/BL1e708LuYxz7iCwG2hSJdFVH5GUU7Br8H7x1kZllpXfjNfOOfffqgysmonkmDoaONrMPjOzL8zszCqLrnqJpKzuAi43s/XAO8ANVRNaXCnv91jMJMU6AKm9zOxyIAs4JdaxVEdmlgD8FRgd41DiQRJek/oQvJad2WbWyzm3K6ZRVU+XAc875x4ys0HAi2bW0zlXGOvApPxUE49cNtAu6HVbf1rIZcwsCa+panuVRFd9RFJOmNlpwG3ACOdcXhXFVt2UVVb1gZ7ALDNbi3dublot7NwWyTG1HpjmnMt3zn0HfIOX1GubSMrqGuA/AM65OUAa3k0/5JCIvseqAyXxyM0FuphZRzNLweu4Nq3EMtOAq/znFwIfOb+XRC1SZjmZWT/gabwEXlvPXUIZZeWc2+2ca+qcy3TOZeL1HxjhnJsXm3BjJpL/val4tXDMrCle8/qaqgyymoikrNYBwwDM7Bi8JL61SqOs/qYBV/q91AcCu51zG2MdVChqTo+Qc+6gmY0BpuP1AH3WObfMzO4B5jnnpgH/xGuaWoXXaeLS2EUcGxGW00SgHvCa3+9vnXNuRMyCjpEIy6rWi7CcpgNnmNlyWUPBnQAAAS9JREFUoAAY65yrba1gkZbVLcAzZvYbvE5uo2tbZcPMXsH70dfU7xtwJ5AM4Jz7O15fgbOBVcB+4OrYRFo2DbsqIiISp9ScLiIiEqeUxEVEROKUkriIiEicUhIXERGJU0riIiIicUpJXEREJE4piYuIiMQpJXERKZWZHevfUznNzOr69+vuGeu4RESDvYhIBMzsXrzhOdOB9c65+2MckoigJC4iEfDH4Z6Ld0/zwc65ghiHJCKoOV1EItMEb7z7+ng1chGpBlQTF5Eymdk04FWgI9DKOTcmxiGJCLqLmYiUwcyuBPKdc/82s0TgczM71Tn3UaxjE6ntVBMXERGJUzonLiIiEqeUxEVEROKUkriIiEicUhIXERGJU0riIiIicUpJXEREJE4piYuIiMSp/w/DUD553mI96QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "-4UT1WedDebi" }, "source": [ "## The Implicit Crank-Nicolson Difference Equation\n", "The implicit Crank-Nicolson difference equation of the Heat Equation is\n", "\\begin{equation}\n", "\\frac{w_{ij+1}-w_{ij}}{k}=\\frac{1}{2}\\left(\\frac{w_{i+1j+1}-2w_{ij+1}+w_{i-1j+1}}{h^2}+\\frac{w_{i+1j}-2w_{ij}+w_{i-1j}}{h^2}\\right)\n", "\\end{equation}\n", "rearranging the equation we get\n", "\\begin{equation}\n", "-rw_{i-1j+1}+(2+2r)w_{ij+1}-rw_{i+1j+1}=rw_{i-1j}+(2-2r)w_{ij}+rw_{i+1j}\n", "\\end{equation}\n", "for $i=1,...9$ where $r=\\frac{k}{h^2}$.\n", "\n", "This gives the formula for the unknown term $w_{ij+1}$ at the $(ij+1)$ mesh points\n", "in terms of $x[i]$ along the jth time row.\n", "\n", "Hence we can calculate the unknown pivotal values of $w$ along the first row of $j=1$ in terms of the known boundary conditions.\n", "\n", "This can be written in matrix form \n", "\\begin{equation} A\\mathbf{w}_{j+1}=B\\mathbf{w}_{j} +\\mathbf{b}_{j}+\\mathbf{b}_{j+1} \\end{equation}\n", "for which $A$ is a $9\\times9$ matrix:\n", "\\begin{equation}\n", "\\left(\\begin{array}{cccc cccc}\n", "2+2r&-r& 0&0&0 &0&0&0\\\\\n", "-r&2+2r&-r&0&0&0 &0&0&0\\\\\n", "0&-r&2+2r &-r&0&0& 0&0&0\\\\\n", "0&0&-r&2+2r &-r&0&0& 0&0\\\\\n", "0&0&0&-r&2+2r &-r&0&0& 0\\\\\n", "0&0&0&0&-r&2+2r &-r&0&0\\\\\n", "0&0&0&0&0&-r&2+2r &-r&0\\\\\n", "0&0&0&0&0&0&-r&2+2r&-r\\\\\n", "0&0&0&0&0&0&0&-r&2+2r\\\\\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c}\n", "w_{1j+1}\\\\\n", "w_{2j+1}\\\\\n", "w_{3j+1}\\\\\n", "w_{4j+1}\\\\\n", "w_{5j+1}\\\\\n", "w_{6j+1}\\\\\n", "w_{7j+1}\\\\\n", "w_{8j+1}\\\\\n", "w_{9j+1}\\\\\n", "\\end{array}\\right)=\n", "\\left(\\begin{array}{cccc cccc}\n", "2-2r&r& 0&0&0 &0&0&0\\\\\n", "r&2-2r&r&0&0&0 &0&0&0\\\\\n", "0&r&2-2r &r&0&0& 0&0&0\\\\\n", "0&0&r&2-2r &r&0&0& 0&0\\\\\n", "0&0&0&r&2-2r &r&0&0& 0\\\\\n", "0&0&0&0&r&2-2r &r&0&0\\\\\n", "0&0&0&0&0&r&2-2r &r&0\\\\\n", "0&0&0&0&0&0&r&2-2r&r\\\\\n", "0&0&0&0&0&0&0&r&2-2r\\\\\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c}\n", "w_{1j}\\\\\n", "w_{2j}\\\\\n", "w_{3j}\\\\\n", "w_{4j}\\\\\n", "w_{5j}\\\\\n", "w_{6j}\\\\\n", "w_{7j}\\\\\n", "w_{8j}\\\\\n", "w_{9j}\\\\\n", "\\end{array}\\right)+\n", "\\left(\\begin{array}{c}\n", "rw_{0j}\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "rw_{10j}\\\\\n", "\\end{array}\\right)+\n", "\\left(\\begin{array}{c}\n", "rw_{0j+1}\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "0\\\\\n", "rw_{10j+1}\\\\\n", "\\end{array}\\right).\n", "\\end{equation}\n", "It is assumed that the boundary values $w_{0j}$ and $w_{10j}$ are known for $j=1,2,...$, and $w_{i0}$ for $i=0,...,10$ is the initial condition.\n", "\n", "The Figure below shows the matrix $A$ and its inverse $A^{-1}$ in color plot form for $r$." ] }, { "cell_type": "code", "metadata": { "id": "KKV6cPkIDebk", "outputId": "d51c5813-a0ff-4a1f-f44d-4d3cac6e27b5", "colab": { "base_uri": "https://localhost:8080/", "height": 297 } }, "source": [ "A=np.zeros((N-1,N-1))\n", "B=np.zeros((N-1,N-1))\n", "for i in range (0,N-1):\n", " A[i,i]=2+2*r\n", " B[i,i]=2-2*r\n", "\n", "for i in range (0,N-2): \n", " A[i+1,i]=-r\n", " A[i,i+1]=-r\n", " B[i+1,i]=r\n", " B[i,i+1]=r\n", " \n", "Ainv=np.linalg.inv(A) \n", "fig = plt.figure(figsize=(12,4));\n", "plt.subplot(121)\n", "plt.imshow(A,interpolation='none');\n", "plt.xticks(np.arange(N-1), np.arange(1,N-0.9,1));\n", "plt.yticks(np.arange(N-1), np.arange(1,N-0.9,1));\n", "clb=plt.colorbar();\n", "clb.set_label('Matrix elements values');\n", "plt.title('Matrix A r=%s'%(np.round(r,3)),fontsize=24)\n", "\n", "plt.subplot(122)\n", "plt.imshow(B,interpolation='none');\n", "plt.xticks(np.arange(N-1), np.arange(1,N-0.9,1));\n", "plt.yticks(np.arange(N-1), np.arange(1,N-0.9,1));\n", "clb=plt.colorbar();\n", "clb.set_label('Matrix elements values');\n", "plt.title(r'Matrix $B$ r=%s'%(np.round(r,3)),fontsize=24)\n", "fig.tight_layout()\n", "plt.show();" ], "execution_count": 4, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvMAAAEYCAYAAADRbnkLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debwcVZ338c+XEBISliC4kEWDggvDOCwRcFeWCaCSZxidAUWE0ckwAwpu4/ogMs8zI26jzPiAUVRAARG3iEh0FGR0ZAkSgQBq2CQBRbawCWT5PX9UNWk6fW9Xd5/qrq77fb9e9bq91D11+t6kvueeOnWOIgIzMzMzMxs9mwy7AmZmZmZm1hs35s3MzMzMRpQb82ZmZmZmI8qNeTMzMzOzEeXGvJmZmZnZiHJj3szMzMxsRLkxP8Ik3SopJL1q2HUxMzNrJumSPKOOHHZdzOqs9o15SV/JTyYhaY2kp3XYf0HT/klPQpLmSjpR0vGpyhwmSZMl3ZX/nNZJmjPsOvVC0iaSXi3pnyWdJ+mWpt//0QmPsVDSLyTdL+lBSVdLeq+kzVIcw8xGU8Vy6sg8p3ZNVeagSZqeZ1K02R6XtErSt0epI0zSlpIOlvQvkn4g6e6mz/T8RMd4hqTPSrpJ0qOS/iDpe5L2TVG+lWfTYVdgwDYF3gh8Zpx93lLi8ecCHwFu61CHom4CHgUeSVBWLw4Cnpo/3gQ4HPi3IdWlH1sBPymrcEmTge+Q/bwAHgfWAbvm2xsk7RMRD5VVBzMbGcPOqSOBVwK3Asv6LOt3wK+B1X2W060XsqGz8g9Nr28CPAWYCfwvYIGkQyPivAHXrxf7At8uq3BJLyTLwW3zlx4AtgNeC7xG0gcj4mNlHd/6U/ue+Sa/y78eMdYOkp4CvAZ4CLh3EJXqR0TsGxHPj4grhlSFRqB8oeX5KHoY+G/g38mC9PcJy/4/ZA35R8mCchowHXgd2b+zFwGfT3g8MxtNtcqpiDgiz6jSGqFj2C3/eldEPKNpexowBVhA9vMTWQfbqLgLuBD4KLAwVaGSNgcWkzXkrwZ2iYitgW2AT5H9nP5V0l+mOqalNZEa878g68neTdKfjbHPocBmwDeBPw2qYqNI0rZkgfIo8F7gZuB5kvYaasV6sxrYKiJeERHviohzgMdSFCzpGcBx+dP3RcQZEbEuMhcAf5e/d1jeM2JmE5dzKo3GEKGrWt/Iz7+Lga/kL+04qEr16XsR8fSIeE1EnAj8KGHZ/wA8i+wPnNdFxHKAiHggIt5DdmVZjOaV9wlhIjXmAc7Kv47V69F4/czxCpH0XEknSPpJPr760Xwc9GWS3p3/ldv6PbcCF+dPn9VmHN+Rzfs2bmyVNEvS/5N0s6THJC1rt1/Ta9s3jaX7v2PU/2VN4wkPGe+zjuMwskC5ICJWA2fnr/fcO5+P04x8/Ogmko6VdEX+s42yxnDmDev1ZZQN/DVZT9BqYFGbY38X+A3ZifKNJdXBzEZH3zklaTtJ/yTpu5JuzO/ReVjS9ZI+LWlmm+85UlKQDbEB+HJLRt3atG/RjGp7A6ykf81fvzvv8GitiyRdlO9zVT5UsRuNnvmNGvNNHsi//rGbgot+9tQiYl1ZZQNvyr+eHRGr2rz/ifzr7pKeV2I9rEcTtTH/JklP+uySngvsBdwOXNKhnLPJLnO9GngG2RCNrfLv/yRwqaQtW77nj8B9+eP1ZOP4mrd2PSzPJRuz+I/A04E1HepFRNxJ9lc2wPskvaT5/bxeZ5L97s+IiG91KnMMjUZ7oxH/tfzroZKm9Fhmg4BvAf8B7A5En+UN06vzr5dGxKNj7PPD/Os+A6iPmVVbipx6P/A54GDgOWT36UwBXgC8E1jW5krgn8iyqJEzD/DkjGrX6O06o3IfIRvOsS3wpTbvHwPMz+t0eEQULRdJmwK75E/bNubzjGp0ZH2zaNktev3slZK3CfbIny4ZY7fL2HDfg2+GraAJ1ZiPiJuBnwOz2PgfZKO342sFemkvB94GzI2IzSNiW2BzshPnb4B5wJNuFImIF7Hh5HF7yzi+Z0TE19sc51PAncBLI2J6RGwBvL7A5/wmcAYwCThL0hZNb38W2IHs5qZ3dCqrHUk7k33G+8nG7xERNwK/JBtj97peym1yCHAA8E9kw1+2ITtZ3txnucOwc/51+Tj7XJ9/fYEklVwfM6uwRDn1O+CDZDeCNjJqCtl5ewnZxAVnN59vIuLrEfEM4H/yl45ryagXtTlOrxm1hqw3+E/AgZL+qfFe3vP78fzp+yLihk7ltXg+MDV//KTGvKQpkvYhu0r+fLI/KHodM9/TZ6+gF5B1oMEYOZX/W/t1/nTndvvYcE2oxnyucWnyzY0X8hPa4S3vjykijomI0yPitqbXHouI75E1QtcCR0qa1mdd1wL7R0Tj5EpErCj4ve8ga7A/m6wBj6T/BRxFdmXgiIh4YMzvHl+jV/6bEdE8tvxrLe/3agvgHRFxakQ8AhARdzXqm1/ebDflWJHtkj7r1q3t8693jLNP470t8s3MJra+cioiTomIf4uIayNibf7auoi4iuzmz+uBPwNe0Wc9e86ovJH+vvzpJyQ9L+9V/ypZ59gPgf/soU67NT2+UtLv8+2PZH88/JhslpZ/IWuI39/DMaDDZ28aNtrLdmKPderF9k2Pi+TU9uPsY0MyERvz55HdtHmIpOn5a68ku/ljaQ+9AE8SEbeQ/XU7jQ034fTqzIj4Q+fd2tbjAbIgWA/8nbL50htjtj8REf/dS7mSJrEhUL7W8va5+fEOUId5kju4h/aXXhseZ+NhSkW3Qc/+0Pg3Nt6Nas1Ti7oxb2al5VTeAdO4efKlfdWyj4zK/SfZlYJpZI34k8iuHtwLHBURvQyxbG7MP71p244NPdAzgIciop8biDt99ofoPacGOU3x9KbHRXLKGVVBE22eeSLifknfA95AdnPimRS88bWZpP3JZiLZk+wv1Y1ueiWby7Yfv+jnmyPiZ5JOBj4AnJq/vAw4oY9i9yf7XKuAn7Yc7w5JF5NdGn4T2TSPvVja6E1qJ+8J2eimKTOzOkiRU8oWEjqWrPd9LlkjrHUY37AzKiQdBVxL1oifl7/1jxExXi/xeBqdaB+PiEbPf2O9j9lk0wP/b+BkSesj4pM9Hmfcz56X22vZZl2ZiD3z0HQJU9nMM68nu3nlnCLfLOkUskuAh5INY9mUrCeh8Vd140aY6W0LKK6ru+zH0FikCrJe8zdHxON9lNcYQvP1McZsNnrrx5wnuYAUn7sqHs6/tvtjr6F5OJYXjjIz6COnJB0KXEN2I+mfk2XRajZkVOO8NPSMyidt+GDTS9/ocxGnRmP+2pbjrImIWyLiI2zIqXf3cZy65NTDTY+L5JQzqoImamP+IrL/iPuQ9VxsCfwgIu7u9I2SDgTeTraC54lkc9ROiYhtGzcKkd0gCxv3gnQrxVRU+5FdmoXs9/2ScfYdl6StyVbNA3hXu7F+bBges2ub2RKKKnMKrkFr9C6N1wPWeO+hiHiw5PqY2WjoKackPZVsIb/JwNfJerunRsQ2TRnVuGo69IzKh24232e1a9PQom7Lmks2CQPAdePs+vP86zMkbdXLsahPTjVfASmSU3eWWBfr0YRszOdDOM4h+/yNudjPGvs7nuQN+dcvRsRHI+KmNuP6np6gmn1TtrBTo3HdOLF9WtJzeizyb9gwS0ARpawIK+klTTc1dbv1OhVnrxoz1Yy1AAxsmB2gr/s1zKw++sipA8mG1FwPvDEirmoztWMlMir3frJOptVkU27uRDZTTC8avfJrGf982jyMs5SpjyW9p4+cek8ZdRrDjWz4GbTNqXyK1Mb88te328eGa0I25nONS5iTyeZ//17B75udf7263ZuSnsXYK8o1hqUMavrBRWRjy5eTje3/Cdll1bPy3pBuNRrnJ5P1foy1HZbv96Z8doLUNuPJNzZ1sz2lhPqMp7FQ2MsljfWH0P751x8PoD5mNjp6yalGRl3TbihkPivOeGtaDCynJO3Ohqkh306WMQH8g6SDeiiycfPripaZ1lo1xubfUeLV0C3oPacGdpNp/vmX5k/3H2O3vYCt88fOqQqasI35fIquE8l6AI7v8B+/WWPhhD8f4/1/ZeyTYGMqyK3HeD8ZZSvuHUI2xvLw/K79I8nq/2Kym2K7KW9HNsx8cG5E3D/WBnyXbBze08kW/kgqIi6JCPW4vSp1fTr4FvAY2ewJb2t9U9LryHo8goL3bJjZxNBjTjUyapcx1q34e7KFpMbSyKkZRevZi/w+gK+S/aFyfkScFREXs2EI0OmStuuy2EZj/tqxdsizrDEj27ldll9YRJzYR06dWFa9xtBY/PFNktpNPdm4UnBVRPy6zfs2ZBO2MQ+QD5N5T0QUnsWGDVN6/YOkv5O0GYCkZ0o6g6xX+r4xvve3ZI3rrSX9dc8V7yAfN/jZ/OlHImIZQETcTjb2EuAESXts/N1jatzQenOjvLHkfzj8IH9aylCb1CRtrWwJ9O3yAGn835je/LrarG6rDUuWX9L6XkT8ng2/i49LenPjqkje8/Tl/L1zIuKa5B/MzEZaDzn1X2SdA7sAp0iaASBpK0nvJVsZ9p5xvr+xcNAh+X1SZTmZbMGi5lXLIbsZdjnZVeXPd1lmY5jNRuPlJW0p6U3ApWRXqG8B/k+X5Q9VS0Zt0/TWjJacal05+Mime9vmtin682QTZWwJXKBsYcjGz+zjbFjw8oNtvtcqYEI35nv0FbKljTcFTgcekXQf2X+EI8guGbZtlEXEw2zofT1f0v2Sbs23JCvH5f+JzwS2IrvJ5+SWOnwV+AZZb8hZee9IpzLFhsVLio45byyRfbCkbcbdsxq+S3azWWObk7/+yZbXD2v73eP7MNlKuZuT/W4elvQw8H2y5cyvBI7up/JmZgB5z+ln8qfHAvflGXUf2cqqPwZOG6eIs8jW8ngZcLekVXlG/SxVHSX9JRs6lv4uIp5Y/yO/+nA4WcfXIflV5iJlbsuG8/a7W8ag30O2YvlXyaaSvgLYJyLG6nirquYs+mXT679oee+Z3RSad8AtIPsjb3dguaTVZD+z95L9cfiBiPhhvx/AyuHGfJfyaR33Az4G3Ew2vnAtWY/96yLiXzoUcTTwb2Q3nUwhm2nmWaQbI/de4OVk00cdMcb0kUeT9Ya8gJbG/hheSTZPMRRvzH+fbHjJFOBvC35PLeU3n72O7Od+GdnPJcjm/H8f8DLPYmNmqUTEu4CFZPd2PQZMyh8fD7yGJ98A2vq9N5KNnb6IbMjOM8gyavZY39ONvHPny2TDUf9fRFzUpg7L2DCW/rNj9Ca3al4saiuePAZ9KlmH2zfI5u1/SUTc2tsnqKeI+BX51Ryyts0Ussb998lWuv3YEKtnHai3BdbMzMzMzGzY3DNvZmZmZjai3Jg3s8qSNEnS1ZIuGHZdzMzMACR9SdJdktouTqbMKZJWSLomn4a1NG7Mm1mVHYcX0zIzs2r5CnDAOO8fSLYA2k5k96+cWmZl3Jg3s0qSNJvsZr0vDrsuZmZmDRFxKXDvOLssAM6MzGVk04e2m8M/iTJW5yxku6dMirlzJpd6jN9cM63U8gFiy/KPYZbao4/ex5rHH066wuP8V0+Pe+5dV2jfq655bDnwaNNLiyJiUctunwH+mWzuY7OB2kxTYirTSz3Gc1/4SKnlw2By0Cy1R3mYx+OxqmfUeGYBtzc9X5m/dmcXZRQ2tMb83DmTuWLJnM479mH+zF0779SntXt2s+6SWTUsveI/k5d5z73ruGJJsemNJ23/20cjYt5Y70t6LXBXRFwl6VWJqmhW2FSms5f2LfUYS5aMu/5eEoPIQbPULo8fJy8zZUZVzdAa82ZWLwGsp92yBj15KdmCYweRzRG9laSvRsThHb7PzMxsI4kzqpNVbFjEDLJ1GlaVdTCPmTezRIJ1sb7Q1rGkiA9ExOyImAscCvzEDXkzM+tduowqYDFwRD6rzd7A6ogoZYgNuGfezBIJYC3FxiOamZkNUsqMknQO8CpgO0kryVYsngwQEacBFwIHASuAR4Cjkhx4DG7Mm1kSQbCuhBWlI+IS4JLkBZuZ2YSRMqMi4rAO7wdwTJKDFeDGvJkls570jXkzM7MU6ppRbsybWRIBrKvpidLMzEZbnTPKjXkzS6auvR5mZjb66ppRbsybWRIBpYyZNzMz61edM6rj1JSSviTpLknXjfG+JJ0iaYWkayTtnr6aZjYK1hfczFJxRplZUXXNqCLzzH8FOGCc9w8Edsq3hcCp/VfLzEZNRPB4wc0soa/gjDKzDuqcUR0b8xFxKXDvOLssAM6MzGXADEnbp6qgmY2GbHW9evZ6WHU5o8ysiDpnVIox87OA25uer8xf22ilK0kLyXpGeOYsD9c3qxexDg27EmatesqoqUwbSOXMbFDqm1FFhtkkExGLImJeRMx76raTBnloMytZAOuj2GZWRc0ZNZkpw66OmSVU54xK0T2+CpjT9Hx2/pqZTTB17fWwkeaMMjOgvhmVomd+MXBEPmPA3sDqiNjo8qWZ1Vu2IIcKbWYD5Iwys1pnVMeeeUnnAK8CtpO0EvgIMBkgIk4DLgQOAlYAjwBHlVVZM6u29TF6J0Ebbc4oMyuqrhnVsTEfEYd1eD+AY5LVyMxG0nrE4/heGBssZ5SZFVHnjPKUMmaWTF17PczMbPTVNaPcmDezJBrjEc3MzKqmzhnlxryZJSLWxUBnuzUzMyuovhnlxryZJZGtrlfPE6WZmY22OmeUG/NmlkxdL2Gamdnoq2tGuTFvZklE1PcSppmZjbY6Z9TQGvO/uWYa82fuWuoxltyxrNTyAebPLP0QrN1nj/IPYtanANbUdNovm3ie+8JHWLKk3AwpOwNhUDlY/ucw61edM8o982aWSH17PczMbNTVN6PcmDezJOp8c5GZmY22OmeUG/Nmlsy6mi7IYWZmo6+uGeXGvJklEYh1Ne31MDOz0VbnjHJj3sySWV/T8YhmZjb66ppRbsybWRLZUtn1PFGamdloq3NGuTFvZkkEYk3Uc9ovMzMbbXXOKDfmzSyJCGo77ZeZmY22OmeUG/NmlohYX9Olss3MbNTVN6M6/okiaY6kiyVdL2m5pOPa7CNJp0haIekaSbuXU10zq6og6/Uospml4Hwys6LqnFFFeubXAu+OiF9K2hK4StKPIuL6pn0OBHbKt72AU/OvZjaB1PXmIqss55OZFVbXjOr4qSLizoj4Zf74QeAGYFbLbguAMyNzGTBD0vbJa2tmlRWI9VFsM0vB+WRmRdU5o7oaMy9pLrAbcHnLW7OA25uer8xfu7Pl+xcCCwGmMq27mppZ5dW118Oqr998yst4IqOeOcu3lJnVTV0zqvDZStIWwDeB4yPigV4OFhGLgEUAW+kp0UsZZlZNdZ72y6otRT7BkzNq3l9MdUaZ1UidM6pQY17SZLIT5dci4lttdlkFzGl6Pjt/zcwmiKC+q+tZdTmfzKyIOmdUkdlsBJwO3BARnx5jt8XAEfmsAXsDqyNio0uYZlZv61ChzSwF55OZdaOuGVWkZ/6lwJuBayUty1/7IPBMgIg4DbgQOAhYATwCHJW+qmZWZRFK1ushaSpwKTCF7Dx1fkR8JEnhVifOJzMrJGVGAUg6APgsMAn4YkR8rOX9ZwJnADPyfd4fERcmq0CTjo35iPgZjP9nSkQEcEyqSpnZaEo4P+9jwD4R8VA+jOJnkn6Qz0ZiBjifzKw7qTJK0iTgc8D+ZDfVXylpccu0uB8GzouIUyXtTNaxMDdJBVr4dn0zSyIg2ep6eQPsofzp5HzzDYlmZtaTlBkF7AmsiIibASSdSzYNbnNjPoCt8sdbA3ekOngrN+bNLBF10+uxnaSlTc8X5TOJbCgt6/m4CtgR+FxEtE45aGZmVlDSjGo35W3rYnQnAj+U9HZgOrBfd/Utzo15M0sioJtpv+6OiHnjlhexDthV0gzg25J2iYjr+qymmZlNQKkzqoDDgK9ExKckvRg4K8+x9X2WuxE35s0sicbqesnLjbhf0sXAAYAb82Zm1rXEGVVkytu3kuUWEfGLfGKH7YC7UlWioZ4TbprZUKxnk0JbJ5KemvfII2lzspuMbiy5+mZmVmOpMgq4EthJ0g6SNgMOJZsGt9nvgH0BJL0AmAr8MeHHeYJ75s0siQhYl67XY3vgjHzc/CZkMwJckKpwMzObWFJmVESslXQssIRs2skvRcRySScBSyNiMfBu4AuS3kk2yufIfHKH5IbWmI8tp7F2zz1KPcb8maUWD8CSO5Z13qlPg/gca/cp93dhE0OqS5gRcQ2wW5LCzHrwm2umMX/mrqUeYzD5Ue5ngPp8Dqu/lENB8znjL2x57YSmx9eTrYVROvfMm1kS2XhEj9wzM7PqqXNGuTFvZsmM4jLYZmY2MdQ1o9yYN7MkArF2feFpv8zMzAamzhnlxryZJZNwdT0zM7Ok6ppRbsybWRKJZ7MxMzNLps4Z5ca8mSVT9ZuLJG0DzMlnyzEzswmk6hkFveWUG/NmlkRZK8D2S9IlwMFk57urgLsk/Twi3jXUipmZ2cBUNaOg/5yq/p8oZjYy1qNC24BtHREPAIcAZ0bEXsB+g66EmZkNV0UzCvrMKTfmzSyJIFuQo8g2YJtK2h74G8CryJqZTUAVzijoM6c6NuYlTZV0haRfSVou6aNt9pki6euSVki6XNLcbitiZqNvfWxSaBuwk8iW3L4pIq6U9Gzgt4OuhJXDGWVmRVU0o6DPnCoyZv4xYJ+IeEjSZOBnkn4QEZc17fNW4L6I2FHSocDJwN8W/wxmNuoixNoK3lwUEd8AvtH0/Gbgr4dXI0vMGWVmHVU1o6D/nOr4qSLzUP50cr5Fy24LgDPyx+cD+0qq5l0GZlaaKl7ClPRcST+WdF3+/IWSPjzQSlhpnFFmVlQVMwr6z6lCf6JImiRpGXAX8KOIuLxll1nA7QARsRZYDWzbppyFkpZKWrpmzcNF62hmI6DC4xG/AHwAWAOQT/d16KArYeUpJaN4rOxqm9kAVTijoM+cKtSYj4h1EbErMBvYU9IuPVSUiFgUEfMiYt7kydN7KcLMKqyiJ8ppEXFFy2trB10JK08pGcWUtJU0s6GraEZBnznV1eChiLgfuBg4oOWtVcAcAEmbAlsD93RTtpmNtsYcvhU8Ud4t6TnkQy8kvR64c9CVsPI5o8xsLBXOKOgzp4rMZvNUSTPyx5sD+wM3tuy2GHhL/vj1wE8ionXMopnVXEXn8D0G+DzwfEmrgOOBfxx0JawczigzK6qiGQV95lSR2Wy2B86QNIms8X9eRFwg6SRgaUQsBk4HzpK0ArgXj0c1m3iCSq6ul88KsJ+k6cAmEfHgsOtkSTmjzKyzimYU9J9THRvz+SD83dq8fkLT40eBN3RzYDOrlwDWrq/etF+STmh5DkBEnDSUCllSzigzK6KqGQX951SRnnkzs44a4xErqHnqrKnAa4EbhlQXMzMbggpnFPSZU27Mm1kyUcETZUR8qvm5pE+SrbRnZmYTSBUzCvrPKTfmzSyZId041K1pZFMYmpnZBDIiGQVd5pQb82aWRFT05iJJ17JhRdBJwFMBj5c3M5tAqppR0H9OuTFvZslU9BLma5serwX+kK8CamZmE0hFMwr6zCk35s0skWrdXCTpKfnD1im+tpJERNw76DqZmdmwVCujIF1O1boxv3afPUo/xvyZpR+CJXcsK/0Yg/gcg/h92PAEsK5a035dRVatdmfvAJ492OqYPdn8mbuWfozB5Ic/h1VfBTMKEuVUrRvzZjZAkY1JrIqI2GHYdTAzs4qoWEZBupxyY97MkqnqTAGStgF2Ipu/F4CIuHR4NTIzs0GrakZBfznlxryZJRFU8+YiSW8DjiOb5msZsDfwC2CfYdbLzMwGp6oZBf3nVOUGD5nZqMpuLiqyDdhxwIuA2yLi1cBuwP2DroSZmQ1TZTMK+swp98ybWTJVG4+YezQiHpWEpCkRcaOk5w27UmZmNlgVzSjoM6fcmDezZCp6CXOlpBnAd4AfSboPuG3IdTIzswGraEZBnznlxryZJRFRyWm/iIi/yh+eKOliYGvgoiFWyczMBqyqGQX955Qb82aWTBUvYUo6BTg3Iv4nIn467PqYmdlwVDGjoP+cKvwniqRJkq6WdEGb96ZI+rqkFZIulzS324qY2eiLUKFtwK4CPizpJkmflDRv0BWwcjmfzKyIimYU9JlT3VxvOA64YYz33grcFxE7Av8OnNxNJcxs9AXFTpKDPlFGxBkRcRDZTAG/Bk6W9NuBVsLK5nwys3FVNaOg/5wq1JiXNBt4DfDFMXZZAJyRPz4f2FdSZe8yMLNyRMFtSHYEng88C7hxeNWwlJxPZlZUxTMKesypoj3znwH+GVg/xvuzgNsBImItsBrYtmglzKwGIt0lTElzJF0s6XpJyyUd12u1JH087+E4CbgWmBcRr+u1PKsc55OZdZYwowAkHSDp1/kQvvePsc/fNOXY2eOU1VdOdbwBVtJrgbsi4ipJrypa8BhlLQQWAkyZOqOfosysitJ1aawF3h0Rv5S0JXCVpB9FxPU9lHUT8OKIuDtZ7awSUuZTXt4TGTWVaf0WZ2ZVkyijJE0CPgfsD6wErpS0uDmjJO0EfAB4aUTcJ+lp4xTZV04V6Zl/KXCwpFuBc4F9JH21ZZ9VwBwASZuSTalzT2tBEbEoIuZFxLzJk6f3Ul8zq7D161Vo6yQi7oyIX+aPHyQbDz2rlzpFxOfdkK+tZPkELRnFlPJqbWZDkSqjgD2BFRFxc0Q8Tnb+WdCyz98Dn4uI+wAi4q6xCus3pzo25iPiAxExOyLmAocCP4mIw1t2Wwy8JX/8+nyfik4AZGZlCLq6hLmdpKVN28Kxys1nH9kNuHwgH8RGhvPJzIpKnFFPDN/LrWTjDqfnAs+V9HNJl0k6oKzP1vM885JOApZGxGLgdOAsSSuAe8lOqmY2kQRQfBaAuyOi49RbkrYAvgkcHxEP9FE7m0CcT2a2kRIyqoNNgZ2AVwGzgUsl/XlE3N9nuW0PVFhEXAJckj8+oen1R4E3pKyYmY2elP2dkiaTNeS/FhHf6qOc5wArI+KxfFY21E0AABx6SURBVFz1C4Ezyzih2vA4n8ysk4QZ9cTwvdzs/LVmK4HLI2INcIuk35A17q9sLazfnKrmurZmNpoSzfuVTx14OnBDRHy6z1p9E1gnaUdgEdkJeMxZBczMrKbSzU15JbCTpB0kbUZ2xW9xyz7fIeuVR9J2ZMNubh6jvL5yyo15M0sk6YIcLwXeTHZD47J8O6jHiq3PpyT8K+A/IuK9wPY9lmVmZiMpXUblmXIssIRsgobzImK5pJMkHZzvtgS4R9L1wMXAeyOi7c339JlTPY+ZNzPbSKJLmBHxMyDVwj5rJB1GdhNkY97eyYnKNjOzUZFwKGhEXAhc2PJa8xC/AN6Vb530lVPumTezNAJivQptA3YU8GLg/0bELZJ2AM4adCXMzGyIqptR0GdOuWfezBIaykmwk/0j4h2NJ/mJ8tFhVsjMzIahkhkFfeaUe+bNLJ10Nxel9JY2rx056EqYmdmQVTOjoM+ccs+8maVToaV48vGHbwR2kNQ8y8CWZPONm5nZRFKhjIJ0OeXGfJ/W7rNH6ceYP7P0Q7DkjmWlH2MQn2MQvw8bQ3cLcgzC/wB3AtsBn2p6/UHgmqHUyGzA5s/ctfRjDCY//DmsT9XLKEiUU27Mm1kyKReN6ldE3AbcRnZTkZmZTXBVyihIl1MeM29m6VRwPKKkQyT9VtJqSQ9IelDSA4OthZmZDV0FMwr6zyn3zJtZMhrOlF6dfBx4XUTcMOyKmJnZ8FQ0o6DPnHJj3szSGN4sAJ38wQ15M7MJrroZBX3mlBvzZpaIqnhzEcBSSV8HvgM81ngxIr41vCqZmdlgVTajoM+ccmPezNKpZq/HVsAjwF82vRaAG/NmZhNJNTMK+swpN+bNLJ0Knigj4qhh18HMzCqgghkF/eeUZ7Mxs3QqOFOApOdK+rGk6/LnL5T04cHWwszMhq6CGQX951ShxrykWyVdK2mZpKVt3pekUyStkHSNpN2LfwQzq4XGghxFtsH6AvABYA1ARFwDHDroSlg5nE9mVkh1Mwr6zKluhtm8OiLuHuO9A4Gd8m0v4NT8q5lNIFo/7Bq0NS0irpCedIJeO6zKWCmcT2bWUUUzCvrMqVTDbBYAZ0bmMmCGpO0TlW1m1o+7JT2H/OKppNeTLZ9tE4Pzycyqrq+cKtozH8APJQXw+YhY1PL+LOD2pucr89ccmGYTiKp5c9ExwCLg+ZJWAbcAhw+3SpaQ88nMCqloRkGfOVW0Mf+yiFgl6WnAjyTdGBGXdltTSQuBhQBTps7o9tvNrOoqOIdvRNwM7CdpOrBJRDw47DpZUknyCZ6cUVOZlrKOZlYFFcwo6D+nCjXmI2JV/vUuSd8G9gSaT5argDlNz2fnr7WWs4jsLw+23Gp2df8+MrPuVXR1PUkzgCOAucCmjTGJEfGOIVbLEkmVT3kZT2TUVnpKBf81m1nPKppR0H9OdRwzL2m6pC0bj8kmtL+uZbfFwBH5rAF7A6sjwpcwzSaaak77dSHZCfJa4KqmzUac88nMulLNjII+c6pIz/zTgW/nfyVsCpwdERdJOhogIk7LK3EQsIJsBSsv0mI2AVV0POLUiHjXsCthpXA+mVlhFc0o6DOnOjbm83E8f9Hm9dOaHgfZ4H0zm8iqOe3XWZL+HrgAeKzxYkTcO7wqWQrOJzPrSjUzCvrMqW7mmTczG5Oisr0ejwOfAD7EhguoATx7aDUyM7OBqnBGQZ855ca8maVTzZkC3g3sOM6iQmZmNhFUM6Ogz5xyY97M0qlmr0djrLSZmU1k1cwo6DOn3Jg3s2QqegnzYWCZpIt58lhET01pZjaBVDSjoM+ccmPezNKp5onyO/lmZmYTWTUzCvrMKTfmzSyNit5cFBFnSNoceGZE/HrY9TEzsyGoaEZB/znVcdEoM7PC1hfcBkjS64BlwEX5810lLR5sLczMbOgqmFHQf065Z34ErN1nj9KPMX9m6YdgyR3LSj/GID7HIH4fo6qivR4nAnsClwBExDJJnpbSLJH5M3ct/RiDyQ9/jrqraEZBnznlnnkzq7s1EbG65bXqLh1iZmYTTV855Z55M0unmr0eyyW9EZgkaSfgHcD/DLlOZmY2aNXMKOgzp9wzb2ZpxIYV9jptA/Z24M/Ipvs6B3gAOH7gtTAzs+GpbkZBnznlnnkzS6eCvR4R8QjZEtkfGnZdzMxsiCqYUdB/Trkxb2bpVOhEKel7jFOjiDh4gNUxM7Nhq1BGQbqccmPezJIQlZsp4JPDroCZmVVDBTMKEuWUG/NmlkaAEs4RI+lLwGuBuyJil66rE/HTprK8aJSZ2USWPqMOAD4LTAK+GBEfG2O/vwbOB14UEUufVKVEOeUbYM0snSi4FfMV4IB+q+RFo8zMDEiWUZImAZ8DDgR2Bg6TtHOb/bYEjgMu71BeXzlVqDEvaYak8yXdKOkGSS9ueV+STpG0QtI1knYvWgEzq5GEjfmIuBS4N0GtTiRbjOP+vNxlwA4JyrUKcD6ZWWHpMmpPYEVE3BwRjwPnAgva7PcvwMnAox3KO5E+cqpoz/xngYsi4vnAXwA3tLx/ILBTvi0ETi1aATOrjy6m/dpO0tKmbWGJ1Wq3GEf1Rk5ar5xPZlZIwoyaBdze9Hxl/tqGY2UdB3Mi4vsFqtZXTnUcMy9pa+AVwJEA+V8gj7fstgA4MyICuCzvKdk+Iu4sWhEzq4HiTeS7I2JeiTVp5kWjasr5ZGZdGVBGSdoE+DT5uamA0heN2gH4I/BlSVdL+qKk6S37dPwLxcxqrujly+EuGnU2sBovGlUXziczKyZtRq0C5jQ9n52/1rAlsAtwiaRbgb2BxZLG+gOhr5wq0pjfFNgdODUidgMeBt5f9ADNJC1sXLJYs+bhXoowswqr4up6EfFIRHwoIl6Ubx+OiE7jF200JMsnaMkoHktVRzOriIQZdSWwk6QdJG0GHAo8ccNqRKyOiO0iYm5EzAUuAw5unc2maf++cqpIY34lsDIiGnfink928mzW6S+URmUXRcS8iJg3eXJr54mZjTqtL7YVKks6B/gF8DxJKyW9tcy620hKlk/QklFMSV5ZMxuuVBkVEWuBY4ElZPfpnBcRyyWdJGngCxJ2HDMfEb+XdLuk5+VzX+4LXN+y22LgWEnnAnsBqz0e0WwCStjrHhGHpSvN6sj5ZGZdSZtRFwIXtrx2whj7virdkTdWdNGotwNfyy8l3AwcJelogIg4jezDHASsAB4BjiqhrmZWZcMZD98TSZvlN0va6HM+mVlnI5RR0F1OFWrM5/Ndtg7aP63p/QCOKVxDM6sd5VvVSLoEODIibs2f7wl8gWwaQxtxziczK6KqGQX951TRnnkzs86q2evxb8BFkk4hm8XkQNw7a2Y28VQzo6DPnHJj3sySGfRMNUVExJJ82MWPgLuB3SLi90OulpmZDVgVMwr6z6miK8CamXVWwXnmJf1v4D/IFhc6kWze39cMthZmZjZ0Fcwo6D+n3DNvZmlE8WknB2xbYM+I+BPwC0kXAV8EiiyxbWZmdVDdjII+c8qNeTNLp4KXMCPi+JbntwH7D6k6ZmY2LBXMKOg/p9yYN7NkqjQeUdJnIuJ4Sd+jzSk8Iga+sIeZmQ1PlTIK0uWUG/Nmlk61TpRn5V8/OdRamJlZNVQroyBRTrkxbwCs3WeP0o8xf2bph2DJHctKP8YgPscgfh9lqFKvR0RcJWkSsDAi3jTs+phZ7+bP3LX0YwwmP/w5hqlKGQXpcsqz2ZhZGkVnCRjgyTQi1gHPylcHNTOziaqCGQVpcso982aWTsV6PXI3Az+XtBh4uPFiRHx6eFUyM7OBq2ZGQZ855ca8mSUhKjvt1035tgmwZf5adU/pZmaWXIUzCvrMKTfmzSwZRSXbyNdHxDeaX5D0hmFVxszMhqOiGQV95pTHzJtZGhUdjwh8oOBrZmZWV9XNKOgzp9wzb2bJVGmmAEkHAgcBsySd0vTWVsDa4dTKzMyGpUoZBelyyo15M0unWifKO4ClwMHAVU2vPwi8cyg1MjOz4alWRkGinHJj3sySqVKvR0T8CviVpLMjYs2w62NmZsNVpYyCdDnVccy8pOdJWta0PSDp+JZ9JOkUSSskXSNp914rZGYjrJrjEedKOl/S9ZJubmwDr4Ul53wys65UM6Ogz5zq2DMfEb8GdgXIV6laBXy7ZbcDgZ3ybS/g1PyrmU0UUdlpv74MfAT4d+DVwFH45v9acD6ZWWHVzSjoM6e6DbR9gZsi4raW1xcAZ0bmMmCGpO27LNvMRpjILmEW2QZs84j4MaCIuC0iTgReM/BaWNmcT2Y2pgpnFPSZU92OmT8UOKfN67OA25uer8xfu7N5J0kLgYUAU6bO6PLQZlZ51ZzD9zFJmwC/lXQsWe/tFkOuk6XXVz7BkzNqKtNKqKKZDVU1Mwr6zKnCPfOSNiO72/YbnfYdS0Qsioh5ETFv8uTpvRZjZhVV0V6P44BpwDuAPYA3A28ZeC2sNCnyCVoyiilpKmdmlVHRjII+c6qbnvkDgV9GxB/avLcKmNP0fHb+mplNFMO7cWhcEXFl/vAhsnGIVj/OJzMbX0UzCvrPqW4a84fR/hImwGLgWEnnkt1YtDoiNrqEaWb1VqWbiyQtHu/9iDh4UHWx0jmfzKyjKmUUpMupQo15SdOB/YF/aHrt6PxApwEXkq1gtQJ4BPd+mU1M1er1eDHZWOlzgMvJ7n+ymnE+mVlh1cooSJRThRrzEfEwsG3La6c1PQ7gmF4qYGY1EaD1lTpTPoOskXcY8Ebg+8A5EbF8qLWypJxPZlZI9TIKEuWU51o2s2SqdHNRRKyLiIsi4i3A3mQ9s5fkMwWYmdkEU6WMgnQ51e3UlGZmY6tYp4ekKWRz9R4GzAVOYeNFhczMbCKoWEZBmpxyY97MkmgsyFEVks4EdiEbM/3RiLhuyFUyM7MhqVpGQbqccmPezNKIqNqCHIcDD5PN3/sO6Yn7ikQ2lHqrYVXMzMwGrHoZBYlyyo15M0umSr0eEeF7gszM7AlVyihIl1NuzJtZOhU7UZqZmT2hphnlxryZpRGgdTU9U5qZ2WircUa5MW8Ds3afPUo/xvyZpR+CJXcsK/0YZX8OxSPlFFzP86SZTQDzZ+5a+jEGkx+j/zn2nO+M6oYb82aWTNXGI5qZmTXUNaN8g5iZpdOYLaDTVoCkAyT9WtIKSe8vueZmZlZ3A8woSe+SdL2kayT9WNKzkn+enBvzZpZMqtX1JE0CPgccCOwMHCZp53Jrb2ZmdTbgjLoamBcRLwTOBz6e9tNs4Ma8maURXWyd7QmsiIibI+Jx4FxgQfpKm5nZhDDgjIqIiyOeuEHtMmB2/x+iPY+ZN7MkstX1Cg9I3E7S0qbniyJiUdPzWcDtTc9XAnv1V0MzM5uohpxRbwV+UPTg3XJj3syS6WLar7sjYl6ZdTEzM2s2jIySdDgwD3hlivLacWPezNIofnmyiFXAnKbns/PXzMzMujeEjJK0H/Ah4JUR8Viyo7fwmHkzS6TgLAHFLnNeCewkaQdJmwGHAotLrb6ZmdXYYDNK0m7A54GDI+Ku5B+nSaHGvKR3Slou6TpJ50ia2vL+FElfz6fnuVzS3DIqa2bVlmqmgIhYCxwLLAFuAM6LiOXl1t5GkfPJzIoqO6MknSTp4Hy3TwBbAN+QtExSaR1SHYfZSJoFvAPYOSL+JOk8sr9AvtK021uB+yJiR0mHAicDf1tCfc2syorfXFSgqLgQuDBZgVY7zicz60rJGRURJzQ93i/ZwTooOsxmU2BzSZsC04A7Wt5fAJyRPz4f2FeS0lTRzEZCgNYX28wScj6ZWWc1zqiOjfmIWAV8EvgdcCewOiJ+2LLbE1P05JceVgPbtpYlaaGkpZKWrlnzcL91N7OqSbi6nlknKfMJWjKK0u5VM7NhqWlGdWzMS9qGrGdjB2AmMD2fZqdrEbEoIuZFxLzJk6f3UoSZVZjWR6HNLIWU+QQtGcWUVNU0s4qoa0YVGWazH3BLRPwxItYA3wJe0rLPE1P05Jc6twbuSVlRMxsBNe31sMpyPplZcTXNqCKN+d8Be0ualo8z3Jfszt1mi4G35I9fD/wkYgR/GmbWuwDWF9zM0nA+mVkxNc6ojrPZRMTlks4HfgmsBa4GFkk6CVgaEYuB04GzJK0A7iWbTcDMJhAR3SyVbdY355OZFVXnjCq0AmxEfAT4SMvLzdPvPAq8IWG9zGwU1fREadXlfDKzwmqaUYUa82ZmhdT0RGlmZjVQ04xyY97M0miMRzQzM6uaGmeUG/NmlozW1/RMaWZmI6+uGeXGvJklMppTepmZ2URQ34xyY97M0ghqe6I0M7MRV+OMcmPezNKp5xVMMzOrg5pmlBvzZpZMXefwNTOz0VfXjHJj3mpl7T57lH6M+TNLPwRL7lhWavl7zn+knIJreqI0M0th/sxdSz9G2fkB5X+O38Q95RRc04xyY97M0ghgfT1PlGZmNuJqnFFuzJtZIvWdKcDMzEZdfTPKjXkzS6emc/iamVkN1DSj3Jg3szRqfAnTzMxGXI0zyo15M0skIOrZ62FmZqOuvhnlxryZpVPT8YhmZlYDNc0oN+bNLI0aX8I0M7MRV+OMcmPezNKpaa+HmZnVQE0zapMiO0k6TtJ1kpZLOr7N+5J0iqQVkq6RtHv6qppZ5UUU28wScT6ZWWE1zaiOPfOSdgH+HtgTeBy4SNIFEbGiabcDgZ3ybS/g1PyrmU0UEbBu3bBrYROI88nMCqtxRhXpmX8BcHlEPBIRa4GfAoe07LMAODMylwEzJG2fuK5mVnU17fWwynI+mVlxNc2oIo3564CXS9pW0jTgIGBOyz6zgNubnq/MX3sSSQslLZW0dM2ah3uts5lVVU1PlFZZyfIJWjKKx0qpsJkNUU0zquMwm4i4QdLJwA+Bh4FlQE/XKSJiEbAIYMutZo/eT8vMxhG1nSnAqillPuXlPZFRW+kp/sdsViv1zahCN8BGxOkRsUdEvAK4D/hNyy6reHJvyOz8NTObKAIi1hfazFJxPplZITXOqKKz2Twt//pMsvGIZ7fsshg4Ip81YG9gdUTcmbSmZlZ966PYZpaI88nMCqtpRhWdZ/6bkrYF1gDHRMT9ko4GiIjTgAvJxiquAB4BjiqjsmZWcSM41tBGnvPJzIqpaUYVasxHxMvbvHZa0+MAjklYLzMbNTWe9suqy/lkZoXUOKO8AqyZJRPrR2+soZmZTQx1zSg35s0skdGc0svMzCaC+maUG/NmlkYwkjcOmZnZBFDjjHJj3szSGcEpvczMbIKoaUYVmprSzKyTAGJ9FNr6JekNkpZLWi9pXv+1NzOzOkudUZIOkPRrSSskvb/N+1MkfT1//3JJc9N+og3cmDezNCKyXo8iW/+uI5tT/NIUhZmZWc0lzChJk4DPAQcCOwOHSdq5Zbe3AvdFxI7AvwMnJ/5ET/AwGzNLJgY07VdE3AAgaSDHMzOz0Zcwo/YEVkTEzQCSzgUWANc37bMAODF/fD7wn5KUT5eblEoos9iBpT8Ct3X5bdsBd5dQnUGV72P4GFU5xrMi4qkpKyDporweRUwFHm16vigiFvVwzEuA90TE0m6/12w8PWRUFf+f+xg+RhWO0Uv5lc4oSa8HDoiIt+XP3wzsFRHHNu1zXb7Pyvz5Tfk+yX/WQ+uZ7+WXJGlpRJQ2Prbs8n0MH2OUj9FJRByQsjxJ/wU8o81bH4qI76Y8llmrbjOqLv/PfQwfY9TKLyp1RlWJh9mYWSVFxH7DroOZmVkbq4A5Tc9n56+122elpE2BrYF7yqiMb4A1MzMzMyvuSmAnSTtI2gw4FFjcss9i4C3549cDPyljvDyMXmO+6zG1FSvfx/AxRvkYlSHpryStBF4MfF/SkmHXySa8uvw/9zF8jFErf+AiYi1wLLAEuAE4LyKWSzpJ0sH5bqcD20paAbwL2Gj6ylSGdgOsmZmZmZn1Z9R65s3MzMzMLOfGvJmZmZnZiKpcY17SlyTdlc/P2e59STolXx73Gkm7d1n+HEkXS7o+Xw7+uBKOMVXSFZJ+lR/jo232SbLMr6RJkq6WdEEZx5B0q6RrJS2TtNFc3v3+rPIyZkg6X9KNkm6Q9OKUx5D0vLz+je0BSceX8Dnemf++r5N0jqSpLe+n+H0cl5e/vPUzpPocZtZe2fmUl1GbjCo7n/JySs0o51NXx3A+DUtEVGoDXgHsDlw3xvsHAT8ABOwNXN5l+dsDu+ePtwR+A+yc+BgCtsgfTwYuB/Zu2eefgNPyx4cCX+/x5/Uu4Gzggjbv9X0M4FZgu3He7+tnlZdxBvC2/PFmwIzUx2gqaxLwe7IFKVL+zmcBtwCb58/PA45M+fsAdgGuA6aRTSv7X8COZf2svHnz9uSt7HzKy6hNRpWdT/n3lppRzqfCx3A+DXGrXM98RFwK3DvOLguAMyNzGTBD0vZdlH9nRPwyf/wg2V3IsxIfIyLiofzp5HxrvdN4AdlJArJlfveVulubXtJs4DXAF8fYpe9jFNDXz0rS1mQBeTpARDweEfenPEaLfYGbIqJ1ZccUx9gU2FzZfLLTgDvaHKOf38cLyE5+j0R2J/1PgUNK+Bxm1kbZ+ZQfoxYZVZF8ahynp5+V88n5NCoq15gvYBZwe9PzlWx8oiskv4y0G1mvRNJj5JcXlwF3AT+KiDGPkf/DXw1s280xgM8A/wysH+P9FMcI4IeSrpK0cLxj5Lr9We0A/BH4cn459ouSpic+RrNDgXPavN7XMSJiFfBJ4HfAncDqiPjhWMfo8fdxHfBySdtKmkbWyzGnZZ+UPysz607S/38jnlGDyCcoN6OcT8U5n4ZoFBvzSUjaAvgmcHxEPJC6/IhYFxG7kq0KtqekXVKWL+m1wF0RcVXKctt4WUTsDhwIHCPpFYnL35TssvWpEbEb8DAlzcWqbGGHg4FvlFD2NmS9DjsAM4Hpkg5PeYyIuAE4GfghcBGwDFiX8hhmVg2jnFEDzCcoN6OcTwU5n4ZrFBvzRZbQHZekyWQnya9FxLfKOEZDfknuYuCAsY6h3pb5fSlwsKRbgXOBfSR9NfExGn/RExF3Ad8G9hzrGLluf1YrgZVNvULnk508Ux6j4UDglxHxhzbv9XuM/YBbIuKPEbEG+BbwkrGO0cfv4/SI2CMiXgHcRzaetu0xcj3/2zWzriX5/1eDjBpIPuX1LzOjnE9dcD4Nzyg25hcDR+R3Re9NdrnozqLfnI8BOx24ISI+XdIxnippRv54c2B/4MY2x+h5md+I+EBEzI6IuWSX5n4SEa1/afd1DEnTJW3ZeAz8JdmltNZj9PyziojfA7dLel7+0r7A9SmP0eQw2l/CTHGM3wF7S5qW/xvbl2ysa+sx+lraWdLT8q/PJBuPeHbiz2Fmvev7/18dMmoQ+ZTXvdSMcj45n0ZGVOAu3OaN7B/zncAasr+K3wocDRydvy/gc8BNwLXAvC7LfxnZGLtryC4DLSMb25XyGC8Ers6PcR1wQv76ScDB+eOpZJfTVgBXAM/u42f2KvLZAlIeA3g28Kt8Ww58KH892c8qL2NXYGn+8/oOsE0Jx5hO1suwddNrqY/xUbJAvA44C5iS+ncO/DdZmPwK2LeMz+HNm7f2W9n5lJdRq4wqK5/yMkrPKOdTV8dwPg1pU/4DNjMzMzOzETOKw2zMzMzMzAw35s3MzMzMRpYb82ZmZmZmI8qNeTMzMzOzEeXGvJmZmZnZiHJj3szMzMxsRLkxb2ZmZmY2ov4/hTmUejCYNcAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "cyrIQpUrDebl" }, "source": [ "## Results\n", "To numerically approximate the solution at $t[1]$ the matrix equation becomes \n", "\\begin{equation} \\mathbf{w}_{1}=A^{-1}(B\\mathbf{w}_{0} +\\mathbf{b}_{1}+\\mathbf{b}_{0})\\end{equation}\n", "where all the right hand side is known. \n", "To approximate solution at time $t[2]$ we use the matrix equation\n", "\\begin{equation} \\mathbf{w}_{2}=A^{-1}(B\\mathbf{w}_{1}+\\mathbf{b}_{2} +\\mathbf{b}_{1}). \\end{equation}\n", "Each set of numerical solutions $w[i,j]$ for all $i$ at the previous time step is used to approximate the solution $w[i,j+1]$. \n", "The left and right plot below show the numerical approximation $w[i,j]$ of the Heat Equation using the BTCS method for $x[i]$ for $i=0,...,10$ and time steps $t[j]$ for $j=1,...,15$." ] }, { "cell_type": "code", "metadata": { "id": "d6oJSy1IDebn", "outputId": "0dbd689a-33a5-4cfb-8022-1246fb46a768", "colab": { "base_uri": "https://localhost:8080/", "height": 487 } }, "source": [ "for j in range (1,time_steps+1):\n", " b[0]=r*w[0,j-1]+r*w[0,j]\n", " b[N-2]=r*w[N,j-1]+r*w[N,j]\n", " v=np.dot(B,w[1:(N),j-1])\n", " w[1:(N),j]=np.dot(Ainv,v+b)\n", "\n", "fig = plt.figure(figsize=(10,6));\n", "plt.imshow(w.transpose(), aspect='auto')\n", "plt.xticks(np.arange(len(x)), x)\n", "plt.yticks(np.arange(len(time)), time)\n", "plt.xlabel('x')\n", "plt.ylabel('time')\n", "clb=plt.colorbar()\n", "clb.set_label('Temperature (w)')\n", "plt.suptitle('Numerical Solution of the Heat Equation r=%s'%(np.round(r,3)),fontsize=24,y=1.08)\n", "fig.tight_layout()\n", "plt.show()" ], "execution_count": 5, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAAHWCAYAAACyk9sKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZxkVX3+8c/T3bMxMOyCAoKKCoMS1EHcEFwwowY0LmwaGaPgEoxRMerPSFCTKBp3MHFEQDAuiGJGgSDKIhFRJoKsomwiCiL7MsAw3d/fH+fcmds1VdXV3dVV91Y/73nVa7rrbt+qulV96txzn6uIwMzMzMysqob6XYCZmZmZWTtusJqZmZlZpbnBamZmZmaV5garmZmZmVWaG6xmZmZmVmlusJqZmZlZpbnBaj0l6TxJIWlZn+vYIddRm1y3fj53ko7K2z6x19vuB0lzJX1I0tWSHur2viJp77zOG7u1TrN+mW2fD9YfbrBOgqQTiz9ckv5vgnm/5jfw4JK0gaS/l3SupD9LWi3pdklXSDpN0j9I2qXfdXYi/7E5StIm/a6lQo4FPgLsBKwB/pRvE5L0yvx87j1z5fWHpBvz59pR3Zy3FyTtll+XZdNYxw6lvwET3e7uYvl948+H7pE0JOkFkv5R0imSbijtL2/t4jYOk/QzSXdLuk/SJZLeK2luN7bRLyP9LqDGni7pVRHx3X4XUjM3AdcA9/S7kKmStCPwP8ATSnc/AMwFdsm3VwK/AnbreYGT98/5/xOBVn9kbye9brf0oqB+krQxsCz/+uopvMdfCRySfz6vS2XZ9O1G2tfPJ+3r03UXsLrN9Np+xjXw50P3LALOmamVS5oDfA94Wb5rNTBK2vd3A14r6YURcf9M1TCT3MM6PR+R5OdwEiLiDRGxU0Sc1u9apkLSCOkD4QnArcDbgM0iYsOIWARsDuwHnAQ82LdCuywijsmv2wf6XUsPPJn0Zf4OfyG1Nl4VEVu3uT253wX2yiz7fJiuB4ALgM8AB5P+jnTLv5Aaqw+RvnRvACwE9gXuBHYHvtTF7fWUG1tTcz6witSTdnCfa7HeejHpdQfYNyL+MyLuKiZGxJ0R8f2IOAR4QV8qtOlakP+vZS+EmVXWPcCiiHh+RLw7Ir4BPNyNFUvaGnhn/vV9EfHViBiN5AfA3+ZpB0natRvb7DU3WKfmVuCY/PNRudetY6UxKzu0mN7yhKDyiTeSFkn6hKTrJD0o6XpJH5E0vzT/iySdlcdXPiDpJ5L2nKC+DSX9P0kXS7onn3TyW0mfl7Rdi2XKdW0i6WhJv5a0qjyWSxOcOKTkAEmnS7pV0sOS/pDrfpekzRvm31bSEZL+J9e4StK9eczOh2dg3NVT8/9/ioiV7WaMiIdaTZP0NKVxzr/Pj/H2/Dq9erIF5ec8JJ3XZp71TopQHpNdmq08nqpx3rYnVeRxU2+SdL6kO/M+c4Ok5XkIRbNlxp14JOm5kn6Qn4sHJf1K0uGSNJnno2EbT5D0pfzeeEjSXXlferOk4YZ5l+Xn47x81/YNz8eyCba1d16+GA7wzw3Ltzxpa6qPXdK+kv47v1dWS7pN0vcl/eVEz00VSHqepG9Kujm/D+6Q9CNJB7V67JKeLunjkv5X0k2l5c5r9rrmZQI4If+6V+Proh6ON5b0mPy++EPeJ6+X9Gmlz82m72V1cJJo4/upYdoWkt6e95VfK41rfEDSVXnbj2myzMB+PpTrzbUdLukXSmM+Q9KMDOXKjcexmVg38GpgHqlRvLzJtv8b+A0g6trRFhG+dXgjjeEJ4JukQ7/35N8PbTLv1/K0E5tMi3zbocV2dijmaTLtvDztXcCv88/3k8aqFOtdked9OzBGGsNyT2n6w8BzW2x7Z+DG0ryP5PUXv9/ZbNlSXe8Frss/PwTcC9zdZL5lTdaxMXB2aVtjeXsPlu5b1rDMqQ2P6478eIv7rgW2ncxzPME+8N683Gpg/hT3o8MaaryLdGJP8fvJwHCb57jxOViW7z+vzTaPatwfgc+RvnwV2/1z/r24fa7d8qVpGwBnldazmjTWrfj9QeAVTZbbO0+/MT+GNfk1Ly8bwGen+Dz/VcO+czfj3ydnAwtL8x+QH/edefpow/NxwATbe06er9jm/Q3L39qtxw7MYd1nTHG7p+H3o6fyvE3wGG/M6z5quvMCRzepf6z0+zeAoSbL3V6a5wHS+6e8ntOBkYZlbi09P6sbXxfgOR0+/h1K29l7Cs/fzsBtpXXcTzpaF8BvgXfT5L1MB59X5X2qybR/L23zEdLnZPkz5zZg14ZlBvbzoVTvV0lDvCKvv9iXduv2e6eD98lbp7me4m/hijbzfCHP84tePb6uPlf9LqBON0oN1vx7sdPfBMxrmHemG6x3kxqsz8v3zwXenD+MAvhQ/mD4N2CTPM/2wIWtdlhSg/GGPP0UYFdywwl4PPBfedqtxTqb1HVffj6Wkv/YADs2mW9Zk+3/IE9bBfx9qW6RPug/3PjBBnwUeAfwxNL25gB7Ab/I6zt9Ms/xBPvAC0qv34nARpNc/jmsa6x+m9yYBjYEPsi6P9j/1Oa1X9Zw/zKm0GDtdH/sYPn/ZN0XlLeQ3wvAk4BzWdeweFLDcnuXpj1M+jDdKk/bBPg867647DLJ5/kJrPuidR7w5Hz/PNIXhofytOOaLFvUdeNkPyNi/OfEUW3mmdZjJ41/Kxo5ryU3vIGNSOOq783TD5rKY2hT940TPbZO5iUduiw+Sw4FNs73LyB9cbglT/9Ak2W/DhwIbF26byHw+tJy722y3ITvkw4e0w6l98vek1x2DnBlXvY64Pn5/iHSGMPbWNcYO6/VdjvYp9bbb0mfpx8gHSEayfcNA88gnUAawBWAmiw7iJ8PRb335breBmyQpz2KdNi+a++bDt8n022wXpXX87E287yt9LjXe62rfut7AXW6sX6DdRHpm2oA72yYd6YbrI9QagiWpn+ltP7jm0zfnnWNosc2TPuXfP/X2zwHZ+Z5jmhR12rgKW2WL+Zb1nD/y0ofPku79HptRvojMNb4XLd7jidYp0qPIUiNoh+QviAspaEh32T5H+fl/pfmvaj/VvpAWdQwrdVzt4w+NVjz81g0wN/SZLkNSL3cAZzUMG3v0ra/3GK7l+XpR07ydSreB9eS/xA1TD+stL/t2KKuG6e4351I5w3WST920pezsbxvb9di2QPzsldM5z3UZL03lvb7xl7KxluxXxzVsI5N8v79IPAXLbbzbNYdYZk7ifr2zNu8ocm0Cd8nHax/h9LrducEj/8jDcv+TV7uYfIXqBa1d73BOsFjmse6hvReTaYP4ufDUaV1H9bBczqVW0f7Gd1rsBa9w+9oM88rSvVNqrOlCjePYZ2GiLgX+ET+9QOSFvZw89+OiGub3P+j0s8fa5wYEb8jfUAAPKVh8iH5/0+12e7X8//7tJh+ZkRc0Wb5Vt6Q/z8rIv5nCsuvJyLuJPUoi9Sz2Y11ButSAMZIPTsvJ2V2ngncoZTN+rLGZSVtxroTsT4WEaNNNnE06Rv/hqyLJqmyvyb1EN0KHNc4MSJWse498qpm4wuz9fbV7L/z/437akt5TFsxFvgzuYZGxwF/IO0br+l03TNkso/9DaS6vxURv2+x7KmkhtEukh49/RLXsxDYaoJbq78vrybt3z+KiF81myEifkY62rMpqRewIxFxAamXcodm4zK7bFPaP/5FDfMX+9l3I+KaxpXl2n8yY9W2EBEPk4bHADy3y6uv3OdDgzuA49tMX826DObJ3u6cYk1TVbQ/2qXTlD8LN5zBWmaEc1in7wuk8aRbkQ67tHpjddvlLe6/Lf//EOsapo3+ROql2bS4Q+lkqm3zr2e0GdxfBA83PfkK+FmL+yfyrGLbk11Q0jOBt5Iapduy7o1b1rU/XvmLyiGSPgS8itQz8gxS7/UQ6Vv53pI+HRHvKS36NFJDI0hJE83WfY/SRSmeCzydNF66yp6e/7+gRQMc1uUOLiRFRl3VMP3OiLi+xbJ/yP9v2mJ6M48nDW+BdMhxPRExlk9seR3rHkM/TOWxF1++DpH02jbrnpP/347u52N+OCKOajdDPllm+yaTivpfKKldpM9m+f/taPhcyY+7eO22BOazvscAf2xX4zS9ICLOm8T8xX7W9L1fmvb8KVfUhqSdgMPz+ncgNVgaT1jqdiO/ip8PZSsjYk2riRFxIbD1FNdtXeYG6zRFxCpJ/0YaoP5eSV+MiF4ERrf6A1R8KPwp9wa2m2dO6b5yL8yjOtj+Bi3u/3MHyzazVf7/psksJOkI0jf04oN3lPGB3huT/ph1vfc7Im4CPptvSHosqffon0h/bN8t6SeRzs6E9IcV4J5oH9x8c8P8VVbU+Ic289xc+rnZY7qvzbJF0sKcNvO0qgk6q6ufz/NUHnvxXt0o3ybS6r3aL0X9G9BZbWvnUUpkOYXUc1d4mHQiVvG5tiXpi2Mvj3h1otjP2jWi2+2vUybpQNJRoWJfGiOdhFZEKm1Ier66/ZxV8fOhbKp/r6roAdJwmwVt5im/32oX2+chAd3xJeD3pG9575lg3qoq7wubRoQmuO3QYj2tvkV3ndKlT48mNVaPIeWjzouIzSKHd5MOjcL6PQldFxE3RcRnSOPvikMvf9tk1nkzXUsfNOvhqoKq1jUdxXv1XR28TzXJXsBeKOr/XIf1n1ha9lBSY7U4MXO7iJgfEVuW3vNFg3DG3/N1IGlL4MukRt23gCWkhJNNS8/ZZ4rZZ6iMqr4Pe/b3qgeK/b5dL3kx7f6IaPdFoJLcYO2CPAboo/nXf5C0xQSLFG+SVm/ijVvcP5PK10l/bB+33+wQYiuvJu3DZ0XEOyLiqiaHnbZqstyMiojfkE6qgnQmbKH4Nr8g/xFppRia0em3/+KQVrs/CjO1TxU1tttnti393IsejfI2Oqmrbr0sxXulH+/TbphO/cUQiI9GxBciotw7Rx4DOdHnb78U+1knDYpGaw9bq5Sz3aDVe/ylpB7Uq4CDI+L/IuKRhnlm6nOyip8PHZP0nJxxPJVbr6+SVwyl2KXNPIvz/1fPcC0zwg3W7jmBFFWyEfD+CeYtgvS3bTF9924V1amIuIF1f0he2uvtAxfl/ydzolHx/F3SbGI+Ce5Zzab1wAP5//K1xi8hjV+FFlfBUrqOfXGSyS873NZE+xO036eKmqbSu1LUuIekVod3X5j/f4B0vfGZdj3rnpNWz3Mx1hg6f547VQSDz1RvVTGec+kMrX+mFfXvLand4ctm2r7nSWO/WzXoZvp1mUixn7Ubo7pXi/vvLv082b8bxfyXRZPQ+nyS4gsb7y8ZtM+HyZjLxCcXtrpt1mR9M6kYr79nmy81xcnSP+5BPV3nBmuX5IHbR+Vf3874MaGNihOmXtE4QdI84B+6WlznTsz/HyFpm1YzKen2FaROyv+/RFKnf4iLscJPbTH9g3Q2xq9jkp6idAm8dvNsxboP4UuL+3NqQfGh8r7caGr0PtIf3Pvp/AS0Yn/aRtJ6Z1QrXdms3dm/9+b/p/KafpfUENicFBXVuO0NSBdbgHR29Iwfgstjt4vejXe2+EP5ZmAbWJuH203TeT47cRKp7p0lvaXdjJKmejLKTPo2qXGyKXBkuxmb1N/yPZ/Ht/5Lm9XN9OsykWI/e5WkJzZOlPQcWjRm85j3G/Ovzf5ubE7ap5spnrOntLgq1KGk3OJWBurzYTIi4rwOh600u+3d43K/SxqTvAlN9gVJ+5JOagvSRTlqxw3W7vo6qVt+Ae2/sZ6S/z9U0htzI7UYk3kG3T9Ts1MfJ/VObQFcKGn/cg+IpMdKOoz0rfmVXd72mfkm4DuS3lE0inMDebGkT0kqb7eIYnm5pA8UDRNJW0r6JCko+44u17k36RKFJytdFnPtt2ilS+W+jjQcYGPSB/UxDct/KN//dOCbkrbNy24o6f+xrnf+4zmNYEI5quwX+dcTJT01r3NOPpv6e6QT0Vq5Mv//hjaxMu22XVwG8OOSDivtz08iXXVoR9KYw3aNiW77N1Kj6DHA6ZKenGuaJ+lQUug4wFci4roub7t4PpfORKRURFzFujGHX5T0sWI/ApC0kaSXSPoa3W+MT1tE3EF6bwK8X9KX874CgKQFkvaU9B+kWLqy4j3/IUmvKPbXfAb894Fnsu7oRqPidVksaY9uPJZJ+hbp78M8UhLL82DtZUtfTmpwtHvPF383/knSfrmBjqRnkeIM57ZY7kekRspTgM+XPlcXSXovcCztPycH8fOhbyRtrHSp3C3y8MGiHbawfH/xPDUsW1za/LzGaRFxK+nkb4BPSPqb0vvjZay7NPE3IuKyrj+wXmgMZvWtbTDviZQuHNBinlczPjz4xCbzzCEdAi/meYR1lw28g1K4b5Nlz6NJeHxp+t5MEB7dbh2kD4+rSrWtIZ2Bu4rxj+uQydTV4bY3YXwo/2h+PtpdmvU7pWlF0HhxYYTjaBHiztQvHPCWhuchSGex3ttw34ONz1HDOkYbai5fJvFrTOLSrHnaHg2v0X2kb9tBupJNcVGIZvvjGxvq/h2pN+ffS/Mc1Wb5DYAfltaxmvGXy3yICS692Ob5XsYkQrgblt23Yd8p0iOK339E6dKsk6lrgu1uwboLioySEj1uLK9vuo+ddJWiLzbsc/eQDh2XL2967lQeQ5uabqTJ+2kq85LSNMq13p/fC+XLFt/QsMxmrAuaL/a14rNzTX7Oiu3u3WSb55eWvaN4XYBndfj4dygtP9GFA26l4cIOpDGE5Uuz3kcHl2bNy27KusteF++r4mpuvyNd6avpPgV8umFfuav0PM+6z4d29c70jfF/49rdlrVZtunjJbUtTm94bh8o/f4LanjBgOLmHtbu+y4TjImLNOB9H+CTpDf+GGmnOpE0frFpmHYvRLoYwdNIwxrOJX2wbEz6Y3AZ6dvyy0mNqm5v+25Sz/QhpMbEnaRD+neQ/tD8A7CiYbEDSL2SV5Ma/gJ+SmostjpENp0av0TqHf0gqUf4d6QPiQW5zotIvXs7R8RX26xjd1KP/C2kEyLuIfUevTYiXh+TPDQWET8HnkfqZbqbFFn3G9LhtpdTOmmjybInkA4L/iLPtx3p5LeOTl6JFP79UtJhqAtIf4A3ID03xwFPjXXRXj0TEd8nHTr+Mul9tkGu7X9Jhyf/MiJa9cZNZ7u3k8bOfpd0EsmWpOdz+y5uYzQi3k56zb9Geq7nkYaT3ER6nxxO/y+K0FJE/AvwF6TPlN+yLorqFtK15/+RlHFcXuZO0rj0/2BdHNKDpKMIe8X4RIFmXkVq6N9Aet8Vr8tUzmKf6MIBW5G+WJTrvwrYjfS+uIX02XErqcd8d9qEzUfEXaQM2+WkM8KHSJ85XyB9Jt3cZtl3k/b5S0hfZIfzz//ALP18GES5bbEvKZf8ItZ1WlxKGm72vKhhOkBBuVVuZmZmfSRpGenQ7fnR+zGQZpXmHlYzMzMzqzQ3WM3MzMys0txgNTMzM7NKc4PVzMzMzCrNJ12ZmZmZWaW5h9XMzMzMKs0NVjMzMzOrNDdYzczMzKzS3GA1MzMzs0pzg9XMzMzMKs0NVjMzMzOrNDdYzczMzKzS3GA1MzMzs0pzg9XMzMzMKs0NVjMzMzOrNDdYzczMzKzS3GA1MzMzs0pzg9XMzMzMKs0NVjMzMzOrNDdYzczMzKzS3GA1MzMzs0pzg9XMzMzMKs0NVjMzMzOrNDdYzczMzKzS3GA1MzMzs0pzg9XMzMzMKm2k3wVY70haCnwOGAaOi4iPN0yfB5wEPAO4AzggIm5st865mhfzWTgzBc8iGqrBd8c61Dgy3O8KJhQj1X4ex+ZUuz6AsRr85Rqb2+8KJhZzo98lTGjB3NX9LqGtVbfex+q7H1SvtveXL1gYd9w52pV1/d9lD58VEUu7srIeqMHb3rpB0jBwLLAPcDNwsaQVEXFVabY3AXdFxI6SDgSOBg5ot975LGQPvWimyu4O9eyzZMqGFizodwkT0kYb9ruEiW2yqN8VTGjNFtV+HldtPa/fJUzoga2r/8Vk1WOq3xh8ZLuH+13ChJ6y/R/7XUJb/3vot3q6vTvuHOUXZz22K+safvRvt+jKinqk+l+lrVueCVwbEddHxGrgm8ArGuZ5BfDV/POpwIukGrT2zMzMZoEAxrr0r27cYJ09tgF+X/r95nxf03kiYg1wD7B5T6ozMzMza8FDAmzSJB0GHAYwnw36XI2ZmdlsEYxG/XpHu8EN1tnjD8B2pd+3zfc1m+dmSSPAxqSTr8aJiOXAcoBF2qz6A7XMzMwGQBoSMDv/7HpIwOxxMfBESY+TNBc4EFjRMM8K4JD882uAcyJidr4zzMzMrDLcwzqAmsVXRcQaSYcDZ+X7HwD+W9J84KiIOA44DfjnfMj/z8CzO9zgDDyK7tHc6ufLaKT6b0XNmdPvEiY0Nqf6z2PVY6NiuNrvZ4CofkgAYzWocWi4+v0RI+pOhNMgqeMJU91Q7U9Om7RSfNVLgcXAQZIWA0TEGRHxJOBTwIURsSNwBCnqCtJJVq8C3g58OyKu73X9ZmZm1lwQjEZ3bnXjBuvgmXJ8VUQ8EBH/CzzUu3LNzMzM2nODdfA4vsrMzGxAjRFdudVN9Qd8WeU41srMzKz3AhitYWOzG9zDOngmE19Fu/iqViJieUQsiYglc6j+ZRzNzMys3txgHTyOrzIzMxtQHhJgA6FJfNXxEXGlpI8AKyNiBfAV4GRJ1wJ3khq1AEi6EVgEzJX0SuAlEXFVrx+HmZmZjRdQyzP8u8EN1sE0Rtqv03AXICKOLE2PPM+4nyXtQxoacB+wGnjvhI1VgYarHThY9foAmFv9jFNGavA8jlT/oNHYSLVzTsdq8FehDjXWISt2aKj6eZ4jFa9RfeiprPYzMnOq/+luk9Iuh7XkTcBdOYf1M8DR+f7bgX0j4qmkIQMn96ZqMzMzs9bcYB0808lhvSQi/pjvvxJYIMlnVZmZmVVAEIx26VY3NTiwYpPULId1j1bz5DGvRQ7r7aV5Xg38MiIensFazczMrFMBo/Vra3aFG6y2Hkm7kIYJvKTFdOewmpmZWc94SMDgmVYOq6RtgdOAN0TEdc02MC6H1SMGzMzMeqI4S7obt7pxg3XwTDmHVdImwOnA+yPipz2r2MzMzDogRrt0qxsPCRgweUzqCcA1gEiN0cYc1pOBKyWtJsVXvTQv/jHgKcB/SQK4FdgjIm5rvUVBxWOjNFL93bwONcacOtRY7X0RICoevTU2XP0/ZHWIjIrh6g80HBqufj/biKpdo6r/dhkY1f7ktEnLsVbLgJ2AhcDWkhZHxJG5sQrweuD0iJgLvBk4PN//HmBeRCwAHp+Xv7OX9ZuZmVlzAYxFd2514wbr4JlOrNWqiFiT758PNcy9MDMzG2CzdUiAG6yDp1ms1Tat5skN1CLWCkl7SLoSuBx4a6kBa2ZmZtYX1R+UZj0VET8HdpG0M/BVSWdGxEPleRxrZWZm1nvpeuv16x3tBvewDp5pxVoVIuJq4H7SSVg0TCvFWs3vYulmZmbWzlioK7e6cYN18Ewn1upxuQGLpO1JJ27d2JuyzczMzJrzkIABk2OtDgfOAoaB45vEWn0FOFnStaQUgAPz4s8D3i/pEVKu8Nsj4vb1t2JmZma9NpuHBLjBOpjGSPt12reBiDiyNL24WMa4nyPiZFJD9rHAVcCOE21IEqp4DmvVc2KBetQ4Uv0axyqecQowNlLtPzZjNfirUI8c1n5XMLGhoeoHwQyp2jX2+t0ciNFZenB8dj7qAZZzWI8lXQxgMXCQpMUNs70JuCsidgQ+AxzdMP3TwJkzXauZmZlZJ9xgHTxTzmEFkPRK4Abgyh7Va2ZmZh3ySVc2KKacwyppQ+B9wId7UKeZmZlNQjGG1RcOsNnuKOAzEXF/u5kkHSZppaSVq8dHtJqZmZl1XQ2G19skTSaH9eaGHNY9gNdI+gSwCTAm6aGIOKa8cEQsB5YDbDy8RbVHxJuZmQ0MMRqzs6/RDdbBszaHldQwPRA4uGGeIof1Z5RyWIE9ixkkHQXc39hYNTMzs/5IsT5usNoAyDmsJwDXkBI3zmmSw3oycKWk1cBqUqIAknYArs7Lbk06+erfJ9xo1SOZhmvw5q5BZFRU/XUGYrj647Ki4p+6dYhjqkUHUx0io2pQ45yh0X6X0Jao/nM4KCr+0WmTlWOtlpGuUnUzcLGkxQ05rK8HTo+It0o6EDgcuCBPuy4idutlzWZmZtaZOp4w1Q11+J5qkzOtWCszMzOrpog0hrUbt7qpX8U2kSnHWuVpj5N0iaTzJe2JmZmZWZ95SICV3QI8NiLukPQM4HuSdomIe8szSToMOAxgvhb2oUwzM7PZaWyWDglwg3XwTDnWKicFPAwQEf8n6TrgScDK8sKOtTIzM+u9dOGA2XlwfHY+6sG2NtZK0lxSrNWKhnmKWCsoxVpJ2jKftIWkxwNPBK7vUd1mZmZmTbmHdcDkWKvDgbOAYeD4JrFWXwFOlnQtcCepUQvwfOAjkh4BxoC3RsSdvX8UZmZmtj5fOMAGyxjpyEE6egA0xFpFnmfczxHxHUm/Bb4ELAL+TdLZEW2uvyqhiuecaqQGu3kNMk4ZqfbrDBA1qHGs4lmxtchhrUWN1R8tNTQ0NvFMfTai6tfYS7P5wgGz81EPsHxI/1jSxQAWAwdJWtww25uAuyJiR+AzwNF52RHga6Se1V2AvYFHelS6mZmZWVNusA6e6eSwvgS4LCJ+BRARd0REtS8zYmZmNouMhrpyq5saHCu1SWqWw7pHq3nymNcih/VJQEg6C9gS+GZEfKJxA+NirYYca2VmZtYLgZwSYEb6AvM84HX5/7+W9KLGmSJieUQsiYglc7Wg1zWamZnZDJN0vKTbJF3RYvrrJF0m6XJJF0r6i5msxw3WwTOZHFbKOayk3tifRMTtEbEKOAN4+oxXbGZmZh0Zi6Gu3DpwIrC0zfQbgL0i4qnAR8n57DPFDdbBM+UcVlIU1lMlbZAbsnsBV/WobjMzM2ujuHBAN24TbiviJ6Toy1bTL4yIu/KvF5E6yGaMx7AOmDwm9QTgGkCkxmhjDuvJwJWSVgOrSYkCAC8jjWW9I/8+HzhiwiYy6NgAACAASURBVI2q4t97ahAZFRWPBoN6REZFxSOjoPo1Vj12C+oRa1WH7qBhVT96a85Qtc/7VY+fw6CyJ0y9CThzJjdQg7eUTUaOtVoG7AQsBLaWtDgijsyNVYDXA6dHxFzgzcDhABHxXxGxfUQsIKUNXB8Rl/b8QZiZmdlM20LSytLtsKmsRNILSA3W93W3vPHcwzp41sZaAUgqYq3Kh/ZfARyVfz4VOEaS8rCAwkGkSCwzMzOriC5eOOD2iFgynRVI2hU4DnhpRNwx0fzT4Qbr4JlOrNXtpXkOYP38VjMzM+uTCCpzaVZJjwW+C/xNRPxmprfnBqutR9IewKqIaBVlUcph3bCXpZmZmVkPSPoG6YqXW0i6GfhnYA5ARPwncCSps+uL6dpDrJluj207brAOnsnEWt3cEGtVOBD4RqsNRMRycnzFxiNbVn/UvpmZ2UAQY/TmpKuIOGiC6W8mnQfTE26wDp61sVakhumBwMEN8xSxVj9jfKwVkoaA/YE9e1axmZmZTSiozpCAXnODdcBMJ9ZK0hzg+8CmwOmSToqIj7XdoICqRzINVbw+gJHq5/TEUCWjVMYZG6lBjRV/qesQGTU2XIMDO0PVr3FkeKzfJUxoiGo/j9X/xBkcNfhLbpMxnVgr4LXA3RGxAfAM4C2Sduhh+WZmZtZGry4cUDX1q9gmsjbWKiJWk6KpGs/2fwXw1fzzqcCLlEZMB7Awj2tdQOp9vbc3ZZuZmVk7gRiL7tzqxg3WwdMs1mqbVvNExBqgiLU6FXgAuAW4Cfj3iGh5WTYzMzOzXvAYVit7JjAKPIY0jvUCST8qLkJQcKyVmZlZf9TxcH43uME6eKYTa3Uw8D8R8Qhwm6SfAkuAcQ3WcbFWcxxrZWZm1gsBjM3SlIDZ+agH29pYK0lzSbFWKxrmKWKtYHys1U3ACwEkLQSeBfy6J1WbmZmZteAe1gGTY60OB84ChoHjm8RafQU4WdK1wJ2kRi3AscAJkq4kpXWcEBGX9f5RmJmZ2frE6CwN03KDdTCNkY4cBGlMKhFxZGl65Hkaf14N3J9/HyP11k5AMFzx4Maq58RCLbJiY6QGNQ5X/4O86jmnVa8PqMexwRrksA4PVT+HdWRotN8ltKUe58R6SIANjJzDeizpYgCLgYMkLW6Y7U3AXRGxI/AZ4Oh8/6EAEfFUYB/gU/nKV2ZmZmZ948bI4JlODuti4ByAiLgNuJt00pWZmZlVwGgeFjDdW924wTp4ppPD+itgP0kjkh5HutrVdpiZmVnfRYixGOrKrW48htXKjgd2BlYCvwMuJI+BLXMOq5mZmfWSG6yDZ8o5rDna6l3FTJIuBH7TuIHxOayPqv6ZBWZmZgNitIa9o90wOx/1YJtyDqukDXL+KpL2AdZExFW9KtzMzMxaS7E+6sqtbtzDOmAmyGEdAw4ive4PNMlh/SvgJEnzgMuBfSfcoEBVj2Sqen3UJDJqqPofcHWItRqreI116LypRfTWcPUPPtUh1mqOqh5rZb3iBusAiogzgDMa7v4w6fD+PqQTsS4G9mvoQf0FKWXgCGBFRPyuB+WamZlZRzRrhwS4wTp7rI27ApBUxF2tbbBGxI15WvW/dpuZmc0y6cIBs7Nfd3Y202enTuKuzMzMzCrHPaw2aeNirYYda2VmZtYro7O0r9EN1tmjk7irjoyLtZrrWCszM7NeCOQhATbwOom7MjMzM6scN1hniXwJ1iLu6mrglCLuStJ+AJJ2l3Qz8FrgS5Ku7F/FZmZm1miMoa7c6sZDAgaMpKXA50gZrMdFxMeLaRFxhqT7gc8CH5Z0TUQcWVr8o8CGwLkR8VcdbhGGqx2KGMPVf2OGs2K7og75nFVPpKnFc1iDI6Iaqv5oqWFVv8ahqtfY430xAkbr8AaYARX/6LTJkDQMHAu8FFgMHCRpccNsNwHLgK83WcUngb+ZyRrNzMzMJss9rINlWlmrEfFjSXv3pFIzMzObtNl60pUbrIOlWdbqHn2qxczMzLoopQTMzoPjbrDapI3PYd2oz9WYmZnNHqO9HjhbEbOzmT64upa12k5ELI+IJRGxZO7Qgm6v3szMzGwc97AOlrVZq6SG6oHAwf0tyczMzLoh8BhWGwARsUbSCcA1pLCNc4qsVWBlRKzIh/OPAeYAr5b04YjYRdJuwE+ADYAhSXcAB0fEWW03KqDqsVE1iIxiuPofQDFU/RrH6vA8Vjw2qur1AcRwxaOOqEms1dB6595WzhyN9ruEtkSvX+fZO4Z1dj7qAZVjrZYBOwELga0lLY6IIyOiuKrVD4ElwMnAGyNil3z/KuAZETFCGkqwGvh5L+s3MzMza8Y9rINlyrFWEfGb0s9/lHQbsCVw98yXbWZmZp0Ym6UnXbnBOli6Emsl6ZnAXOC6LtVlZmZm0zSbr3TlBquNI+nRpOECh0RE0wFOjrUyMzOzXnKDdbBMK9ZK0iLgdOCDEXFRq/kiYjmwHGDjeVtV/8wCMzOzATFbT7pyg3WwTDnWStJc4DTgpIg4deZKNDMzs6lIV7qanUMCZmczfUBFxBrgcOAs4GrglCLWStJ+AJJ2l3Qz8FrgS5KuzIvvDzwfWCbp0nzbrQ8Pw8zMzGwc97AOmIg4IycAfA54o6TRiDiyNP1iSQcDnwV2Bd6S7/+apAtIvaxDpJzWZwGXtt+iQBX/tlf1nFggRupQY8VfZ6AOR8qqXmMdOm+q/hwCqAY1Dqn6I7qcw7o+pwTYQMhZrMcC+5BSAi6WtCIirirNdhMpr/WIhsVvAZ4dEQ9L2hC4Ii/7xx6UbmZmZm34Slc2SKaTxbq69Os8PGTEzMzMKsAN1sEzrSxWSduRkgJ2BN7brHfVsVZmZmb9MVtTAmbno7aWIuL3EbErqcF6iKStmsyzPCKWRMSSucMb9L5IMzOz2ShSSkA3bnXjBuvgmVYWayH3rF4B7NmluszMzMymxA3WwbM2izVnqx4IrOhkQUnbSlqQf94UeB5wzYxVamZmZh0LUkpAN24TkXS8pNskXdFiuiR9XtK1ki6T9PRuP94yj2EdMBGxRtIJpIamgHOKLFZgZUSsyGNQjyFFV71a0ocjYhdgZ+BTkgQ8GfhpRFzedoMChqr9vScqXh9ADFX/8EwtahyuQ439rqC9WgyPq0GNqkFk1IiaXn27UoYrXmM/PnF6eDj/RFJb4aQW018KPDHf9gD+g0mcMzNZNXjb22TkWKtlwE7AQmBrSYsj4siIKHpafwgsAU4G3pgbq0TE2Xn86jnAqcCVjes3MzOzwRcRPwHubDPLK0hXx4x8OfdNJD16pupxg3XwrI21yjFVRazVWhFxY0RcBqz31VXSM4CtSI1aMzMzq4gih7UiJ101SyXaphsrbsZDAgbPlGOtJA0BnwJeD7y4+6WZmZnZdHRxSMAWklaWfl8eEcu7tfJuc4PVyt4OnBERN6vN5VbH5bCOOIfVzMyshm6PiCXTWL4rqUSdcoN18ExnB3o2sKektwMbAnMl3R8R7y/PlL+BLQfYeP7W1T+zwMzMbAAElcpQXQEcnq+ouQdwT0TcMlMbc4N18KyNtSI1VA8EDu5kwYh4XfGzpGXAksbGqpmZmfVPJ5FU3SDpG8DepKEDNwP/TEoXIiL+EzgDeBlwLbAKeONM1uMG6+B5MTCXFGt1H/DpcqwVcDepd/TJwMPAvkWslaRRoIixCuCnnWwwhit+7l4doo7qEBlV8ZcZ6lJjtV/rqsduAcRQ9Q/saKjacUwAwzWocY5G+11CW6L6++JURcRBE0wP4O96VI4brIMkR1odCzyfdLLVxcBpABFxZJ5nB2B/4AhgRUScWlrFgxGxWw9LNjMzs05FT3NYK8UN1sGyNtIKII8reQVwVTFDRNyYp1X/q7WZmZmtVcRazUY1OIBmkzDdTLT5klZKukjSK7tbmpmZmdnUuIfVyraPiD9IejxwjqTLI+K6xpnGx1ot6nWNZmZms9Zs7WF1g3WwTCsTLSL+kP+/XtJ5wNOA9RqsjrUyMzPrvYrFWvWUhwQMlrWRVpLmkiKtVnSyoKRNJc3LP28BPJfS2FczMzOzfnGDdYBExBrgcOAs4GrglCLSStJ+AJJ2z3lqrwW+JOnKvPjOwEpJvwLOBT4eEW6wmpmZVUiEunKrGw8JGDxjpBMJAxiFdZFW2QLgNmBr4JAi1ioiLpT0cuA40rCC90j6cZEq0JIEQ9X+3hMVrw+oR1ZsHWqswUtd+ZzTGjyHVP05BIZqkBU7UoMc1qGqB9r04WOxVxcOqBo3WAdIKYd1H3IOq6QVDT2lNwHLSDmsjU4C/jUizpa0Ianxa2ZmZhUQzmG1ATHlHFZJi4GRiDg7z3d/j2o2MzMza6sOB3+sc9PJYX0ScLek70q6RNInc4+tmZmZVYTHsNpsNwLsSYqyugn4FmnowFcaZ3QOq5mZWT841soGw3RyWG8GLo2I63PawPeApzebMSKWR8SSiFgyd2SDaRVsZmZmNhE3WAfLlHNY87KbSNoy//5CnMNqZmZWKR4SYLUXEWsknQBcQwrbOKfIYQVWRsSKfDj/GGAO8GpJH46IXYDnAwuBmyQBzAVO72jDwxX/3lOLOKYa1FiDEc1jtXit+11Be1GDOKY61FiLWCuN9ruECc2peI2it69zMHtTAir+0WmTkU+SWgbsRGp8bi1pcUQcGRFFT+sPgSXAycAbc2OViDg3InaMiAWkE7XuAc7o9WMwMzMza+Qe1sEy5VirBq8BzoyIVTNXqpmZmU1KpCzW2cg9rINlOrFWZQcC3+hKRWZmZtY1Y6grt7pxg9XGkfRo4KnAWW3mOUzSSkkrV4+6E9bMzMxmlocEDJbpxFoV9gdOi4hHWs0QEcuB5QAbL3j0LD04YWZm1lsBtTzDvxvcYB0sa2OtSA3VA4GDJ7mOg4APdLswMzMzmy5fOMAGQA78P5x0OP9q4JQi1krSfgCSdpd0M/Ba4EuSriyWl7QDqYf2/F7XbmZmZtaKe1gHTESckRMAPge8UdJoRBxZmn6xpIOBzwK7Am8pTbtR0n8Bl0saAs4G3hnR5pxEiRiq9veeWmScusauqHrGKVS/xqrXB1CH80XqkMM6pOrXONzjnNPJ6seu6JQAGwg5i/VY4KXAYuAgSYsbZruJlNf69YZlnwM8l9SQfQqwO7DXDJdsZmZmHfKVrmxQTCeLNYD5pKtciXQ1rD/NfMlmZmZmrbnBOniaZbHu0cmCEfEzSecCt5AarMdExNXdL9HMzMwmK8IpAWZI2hHYmRSHBXC2pD0j4oKG+Q4DDgOYP2dRb4s0MzObxZwSYINiOlmsfw1cFBH3R8T9wJnAsxtniojlEbEkIpbMHVk47YLNzMzM2nGDdfCszWKVNJeUxbqiw2VvAvaSNCJpDumEKw8JMDMzq4g0LGD6t7rxkIDB82LSSVPXAPcBny6yWIGVwN2kq1Q9GXgY2FfShyNiF+AZwNOB+4HbSTmu359wixX/2hPDFS8QiOHqH+KpQ9yRa5y+qtcHQB0io4Yaz2mtnpEa1DhHa/pdQlvqQ+zWbB3DWoePJutQKdLq+cCGpBOuTgOIiCMjYgWpF3V/4GTg9RGxeUTsIunlwNOARwGbk068OqrnD8LMzMysgRusg2VtpFVErAaKSKu1IuLGiLgMaPxqvRj4SUSsiYgHgMuApb0o2szMzCYWdCeDtY69tG6wDpZmkVbbdLjsr4ClkjaQtAXwAsafvGVmZmZ9Fl261Y3HsBoAEfFDSbsDFwJ/Bn4GjDab17FWZmZm1kvuYR0s04m0IiL+NSJ2i4h9SBcO+E2L+RxrZWZm1msxey/N6gbrYJlypJWkYUmb5593BXYFfjhjlZqZmdnkzdIxAR4SMEAiYo2kw4GzgGHg+HKkVUSsyIf9TwM2ZXyk1RzgAkkA95ISBKqdJ2JmZmazghusAyYizgDOaLjvSABJS4HPAQ8BH42Ij5fmeUjSccCbgQ2AT0n624j43YQbHap4R/1Q9Q99hGvsjorvikD1a6x6fVCLHNZhVb/GEdUhh7XpqRSV0Y9PxToezu+GOnw0WReUMlpfSoqwOkjS4obZLgGWRMSuwKnAJ3pbpZmZmbVT9ytdSVqY2yST4gbr7NFJRuu5EbEq/3oR6aQtMzMzsymRNCTpYEmnS7oN+DVwi6SrJH1S0o6drMcN1tljshmtbwLObDZB0mGSVkpauXrNqmazmJmZWZcFtUwJOBd4AvABYOuI2C4iHgU8j9Q5drSk10+0Eo9htfXkHWcJsFez6RGxHFgOsPEGj6n+QC0zM7NBEED9xrC+OCIeabwzIu4EvgN8R9KciVbiHtbZo6OMVkkvBj4I7BcRD/eoNjMzMxtARWNV0kcl7SNpvQD3Zg3aRm6wzh4TZrRKehrwJVJj9bY+1GhmZmZt1Pikq+uBg4CVkn4h6VOSXjHRQgUPCRgwpeiqYeC4IrqqlNH6c2AT4E7g85KuAH4cESuAs4HNSDvTDcBNEbFfu+2FIIarfXiiDnFMUYOvjq6xO6peYy2ONlb8OQQYGqp+ZNRQDWKthitfYx9afjUdiBcRJwAnSNoa2B84gnSZ9406Wd4N1gFSiq7ah3RS1cWSVkTEVZAyWiUdAPw8IlZJehuwd26sAhxAymB9S0T8VR8egpmZmbVUz8uqAuSs98XAn4ALgNcAv+x0+Rp8T7VJmFZ0VUT8GLivV8WamZlZdUlaKukaSddKen+T6Y+VdK6kSyRdJullbVa3Oeno792ko7y3T+aKmu5hHSzNoqv2aDN/y+gqMzMzq6AeDQmY6Kht9k/AKRHxH/liRGcAOzRbX0T8dV7vzsBfAudKGo6IjjLf3WCdpSaKrppg2cNI406YP3fjLldmZmZmTUVPL8269qgtgKTiqG25wRrAovzzxsAfW61M0l8BewLPJ51Lcw5paEBH3GAdLJONrtprKtFV5RzWRQudw2pmZlZDW0haWfp9ef77XujkqO1RwA8lvQNYCLy4zfaWkhqon4uIlg3bVtxgHSxro6tIDdUDgYPLM5Siq5Y6usrMzKxmutdNdHtELJnmOg4CToyIT0l6NnCypKdExNp4B0mK5PBWKynmabchN1gHSI6uOgG4BhBwTkRcKekjwMqcBvAtYHvgMklrgIsiYh9JuwE/IaUEDEm6Azg4Is5qv1XBULXP3at67BY4eqtbalFjxV/qOjyHqPoHdmrwlmZO5SOjYIhq19ifl7lnW+3kqO2bSD2nRMTPJM0HtgDKHWLnSvoO8N8RcVNxZ86Efx5wCOnyrSe2K6YOH03WoTxAehmwE6lrfmtJiyPiyFJ01VuARRGxAPhH0pl6AKuAZ0TECGmnXE3KbDUzM7PZZ8ILDgE3AS+CtSdTzQf+3DDPUmAU+IakP0q6StL1wG9JPbSfjYgTJyrGPayDZcIB0hFxbmn+i4DX5/t/U5rnj5JuA7YkxU+YmZlZFfToAEPpgkNnkeKojm9y1PY9wJclvStXtqzx0H5EPAR8EfiipDmkHtgHI2JS7Qs3WAdLV2KtJD0TmAtc19XqzMzMbHp6OCImIs4gRVWV7zuy9PNVwHMnsb5HgFumUosbrLNUq1grSY8GTgYOKQ+abpjHsVZmZmbWM26wDpZpxVpJWgScDnwwIi5qtZHxsVbbVP/sBzMzs0EQVP/MzRnik64Gy4QDpEuxVvuVY63y/KcBJ0XEqT2s2czMzDoU0Z1bP0jaPneaIWmBpI06XdYN1gGSr8lbDJC+mnS5tCslfUTSfnm2TwIbAt+WdKmkokG7P+nqE8vy/ZfmqCszMzOzaZF0KHAqqdMM0lHg73W6vIcEDJiIOEPSGPA54I2SRhsGSL9Y0ruBNwNrgIWSto+Ir0m6gNTLOgTMAZ4FXNp2g6p+zmktMk6H+13BxOqQz1mL17rqz+NQ9Uf5qAY1Dg9VOz8UYGRotN8lTGiuql2j+pEJXP3dv5W/I6UZ/RwgIn4r6VGdLuwG64DJWazHAvuQUgIulrQin8lXuARYEhGrJL0N+ARwAOnMvWdHxMOSNgSuyMtO+hJqZmZmNgPqO4b14YhYLaX6JY0wieZ31b/r2+StzWKNiNVAkcW6VkScGxGr8q8XkbrliYjVpZOw5uH9w8zMzLrjfEn/D1ggaR/g28D3O13YDZLB0yyLdZs284/LYpW0naTL8jqOdu+qmZlZdSi6c+uD95GugnU56aqbZwD/1OnCHhIwizXLYo2I3wO7SnoM8D1Jp0bEnxqWcw6rmZlZrwW1HMOahyteGRE7AV+eyjrcwzp4JpvFul85i7WQe1avAPZsMm15RCyJiCVz5izsWuFmZmY2eCJiFLhG0mOnug73sA6etVmspIbqgcDB5RlKWaxLG7JYtwXuiIgHJW0KPA/4TM8qNzMzszZU55OuNgWulPQL4IHizojYr/Ui67jBOmAiYo2kE4BrAAHnFFmswMqIWAF8C9geuEzSGuCiiNgH2Bn4lNIpfE8GfhoRl0+4zYpHCVU+RojqP4dQlxr7XUEHKl5jHZ5Dx1p1x3CfBjJOxpCq/zz2XPVftlY+NJ2F3WAdMHmcyDJgJ9bFWi0uZ7GSBjv/vBRrtTdARJxNGr/6OdKg6Dt7WbuZmZkNpog4fzrL1+C7tE3SlGOtACQ9A9gK+GGP6jUzM7NORZduPSbpPkn35ttDkkYl3dvp8hM2WCU9SdKPJV2Rf99VUscxBNZzU461kjQEfAo4YsaqMzMzs6mraYM1IjaKiEURsQhYALwa+GKny3fSw/pl4APAI3mDl5FO5LGaK8VafTLf9XbgjIi4eYLlDpO0UtLKR9Y80G5WMzMzs3Ei+R7wl50u08kY1g0i4hfFpbSyNZMtznpmsrFWe5VirZ4N7Cnp7cCGwFxJ90fE+8vLRsRyYDnAog23qe/wbzMzszoJapsSIOlVpV+HSB1mD3W6fCcN1tslPYHcgSzpNaRrzls1TTnWKiJeV5pnGbCksbFqZmZm/VODcIdW9i39vAa4kYZzbNrppMH6d6TetJ0k/QG4AXj9JAq0HsqxVocDZwHDwPFNYq0+SepB/XbuOb+p0xw0MzMzsyk4LiJ+Wr5D0nOB21rMP86EDdaIuB54saSFwFBE3DelMq0nJC0lhf0PAV+OiI8DFLFWkt4NPAb4E3AH8LcR8bs8bZQUZwWdNmIFqNqHJ2K42vVBPbIvXWN3VL3GqtcHMFSDGuuQcTqi0X6XMKG5Fa9xqC9nL/V+k13yBeDpHdzX1IQNVkmbAG8AdgBGirGsEfH3k6nSZl7OYD0W2Id1GawrIuKq0myXkA71FxmsnwAOyNMejIjdelq0mZmZDSxJzwaeA2yZO80Ki0hHgjvSyZCAM0hZnZcDvuREta3NYAWQVGSwrm2wRsS5pfkvwsM7zMzMbObMJQ1DHAE2Kt1/L/CaTlfSSYN1fkS8e+LZrAKaZbDu0Wb+tRms2XxJK0mDoT+eIyfMzMysImow2mScfIWr8yWdWAxBnIpOGqwnSzoU+AFQxB8REb5sZ42VMlj3Kt29fUT8QdLjgXMkXR4R1zVZ9jDgMID58zbuSb1mZmZGbWOtgFWSPgnsAswv7oyIF3aycCdD11eTzir/GfB/+bZy8nVaD0w2g3W/UgYrEfGH/P/1wHnA05ptJCKWR8SSiFgyZ87C7lVvZmZmrXXrKlf96aX9L+DXwOOAD5NirS7udOFOGqzvAXaMiB0i4nH59vipVGozbm0Gq6S5pAzWFeUZShms+5UzWCVtKmle/nkL4LmUxr6amZmZTcPmEfEV4JGIOD8i/hboqHcVOhsScC2waqrVWe/kDNYTgGtIgVPnNMlg/RawPXCZpDXARRGxD7AzcLykbYA5wJ10+LrHULUPT1S9PqhHlFAdauzoK3ifVf55HKr+ADnVYBDfUA1qnFPxyCiAIZ/rvb7q71qtPJL/v0XSy4E/Apt1unAnDdYHgEslncv4MayOtaqYHGu1DNiJdbFWi4sM1uwtwM9LsVZ7A0TEhZJuBd4REWdL2hCnQpiZmVVKDb4LtfIvkjYmHbn/AinW6l2dLtxJg/V7+WbVN+VYK0mLgZGIODvPd3+vijYzM7PBlTvUnhgRPwDuAV4w2XV0cqWrr06hNuuP6cRaPQm4W9J3SQOifwS8PyKqf8zIzMxstqhhD2tEjEo6iHQlzilp2WCVdEpE7C/pctZ/eiIi/mKqG7X+axJrNQLsSUoGuIk01nUZ8JUmy66NtZrnWCszM7PeqWGDNfuppGNI7YsHijsj4pedLNyuh/Wd+f+rgfeW7hfpcp5WPZONtdqrFGt1M3BpaTjB94Bn0aTBGhHLgeUAizbapr5vHTMzM+uV4tLvHyndF3SYFNCywRoRt+Qfd2y8MoGknSZTofXM2lgrUkP1QODg8gylWKul5VirvOwmkraMiD+TdiDn7ZqZmVWEor4nXUXEpMetlrUMWJH0tjwc4MmSLivdbgAum85GbWZExBrgcOAsUs/4KUWslaT98myfJF3T99uSLpW0Ii87ChwB/Di/7gK+3PMHYWZmZq2FunPrMUlbSfqKpDPz74slvanT5dsNCfg66YScjwHvL91/ny/LWmljrLuOxShAQ6zVGcBjSF9W/gy8A0DSC0iN2SLK6kA6SYiQiOGK55xWPfcSZ8V2Sz1qrHj3SA2eQw1VP3FvuAY11iErdm7Fs2JV4wGlfXAicAJpSCLAb0jjWdcbethMuyEB95CiBw6aXn3WKzk24lhgH9blsK6IiPIVqy4BlpRyWD8BHJDjrnbL69mMdMGIH/b0AZiZmVl79W0jbxERp0j6AKy92FHH30hq8F3aJmFtDmtErAaKHNa1IuLciCiuYHUR6cSsRq8BzizNZ2ZmZhVQjGOd7q0PHpC0ObnJLelZpI7RjnRy4QCrj+nksJYdCHy61ULjY602mXyVZmZmNtu8G1gBPEHST4EtSR1kHXGDdZZqksNa3P9o4KmkE7eaGhdrtWjb+h6cMDMzZQvCYgAAIABJREFUq5ua/tWNiF9K2gt4MunE7msi4pFOl3eDdbBMJ4e1sD9w2mR2IjMzM+uBGsdaSZoPvB14HqnZfYGk/4yIhzpZ3mNYB8vaHFZJc0mH9leUZyjlsO7XkMNaOAj4xoxXamZmZrPJScAuwBeAY/LPJ3e6sHtYB0g+4+5w4ALgUcB9wL6SDgBWRsQKUnTVo4DLJQVwPyk14HeSdiDtQMdIGgLOBt4ZES2/zwXVj2QKVbs+qEscU78rmFgdaqx6N0HlY7eAGrylGVH1Y63mVDwyCmCIaj+PfdkVq/8WbeUpEbG49Pu5kq5qOXeDin902hScBTwM7AQ8mtRj+s3cWCUiXgzsB2wWEQuAI1l3qd3HkC4KsSvwFGB3Gsa4mpmZWR9Fl26998ucDACApD2YxBU13WAdPNOJtgpgPjAXmAfMAf7Uk6rNzMxskD0DuFDSjZJuBH4G7C7pckkTXkHVQwIGz5SjrSLiZ5LOBW4hHek4JiKunqlCzczMbHLqetIVsHQ6C7uHdRYrRVt9Mv++I7Azqcd1G+CFkvZsstxhklZKWvnIIw/0smQzMzPrEUlLJV0j6VpJ728xz/6SrpJ0paSvt1pXRPwOuBfYGNi8uEXE7/K0ttzDOnimE23118BFEXF/nudM4Nmkk7jWKuewbuQcVjMzs4HTyeXeJT0R+ADw3Ii4S9Kj2qzvo8Ay4DrWjaIN4IWd1OMG6+BZG21FaqgeCBxcnqEUbbW0IdrqJuBQSR8jDQnYC/hsT6o2MzOzifWum2jtOTEAkopzYspn9h8KHBsRdwG0iMss7A88IZ9fM2lusA6eF5NOmrqGFGv16Yi4UtJHSGfj7Qh8hHRS1S8l3QRcFxH7kQZEP50UdXU7cEpEfH+iDVY9SiiG+13BxKr+HIJr7JbK11iDyKihWkRvVb/GoRrUOFz5eLAeP4e9vXBAJ+fEPAkgX2p1GDgqIv6nxfquADYB2jVqW3KDdYCUuu+fT+6+B04DiIgj8zwvAB4VEaskvQ3YOyIOkPRy4GmkjNZ5wHnAUb1+DGZmZtYTW0gqx0otz0P+JmMEeCKwN2kI4k8kPTUi7m4y78eASyRdQYrfBCB3mHW0IRscE3bfR8S5pfkvAl6ff14M/CQi1gBrcsTEUuCUXhRuZmZmHeheD+vtEbGkzfROzom5Gfh5vpz7DZJ+Q2rAXtxkfV8FjgYuh8lfEaLqB6dscpp132/TZv61kVbAr4ClkjaQtAXwAsbvqGZmZtZvvbtwwISXewe+R+pdJbcdngRc32J9qyLi8zkL/vzi1lEluId11ipFWu0FEBE/lLQ7cCHwZ1Kgb9Pr9kk6DDgMYN78TXpSr5mZmfVO6XLvZ5HGpx5fPicmX0HzLOAl+RKro8B7I+KOFqu8IJ/UvYLxQwJ+2Uk9brAOlulEWhER/wr8a57n68Bvmm3EsVZmZma9J3p74YCIOAM4o+G+/9/evcfLVdb3Hv989849XBIIFwXKRaEagSIGSDneqgJB24hFKaAQ1GMOBapV8IinPSjQU5F6aOlLXmqqIHBERU7VWIKRBqmtniBBAUmQm4AiAYRwKYHcf+eP55lkZZi9Z+/suaw1832/XvPK7DXPWs9v1syePHvNs77rvML9AD6Wb828Nv87u7DMsVZ9apsjrfIJW9Mi4ilJBwMHAz/oWOVmZmbWXEUPE0XEH41lfQ9Ye8gID9//HbAd8C1JAL/OZ+iNJx2uh3QlivflE7DMzMzMxkTSbsDfAi+PiGMlzQT+MCK+MpL1PWDtMcMdvpc0hzRl4D+Bb0TERYU2ayT9I3AmMAn4vKT5xStaNCSIwXIHN8ZAuesDKpF9WYX9WPqMUyDKvhsrkM05MFD2bE4YV4Eax6vhaQqlUvYa1enf587msLbaV4ErSFMSIU07/CYwogFrBT7erRUKGa3HkiKsTsp/3RRdExEHRcQhwMXAJR0u08zMzIbTuZSAlpBUOzg6IyKuJUda5W9xR/wXiQes/WNzRmu+LFoto3WziHiu8ONUKjtTxszMrEdVbMAK/DT/u1rSzrXeJc0Gnh3pRjwloH+M5BJrSDqTdLbfBEZ45p6ZmZnZEGoTJz5GirR6Rb6U6y7Au0e6EQ9YbSsRcRlwmaSTgb8G5tW3cQ6rmZlZd1RwDusukmqxV98mnWcjUhbr24A7R7IRD1j7x4gyWgu+AXyh0QNb5bDu6BxWMzOzjqne/7qDpHSi+lPUpoxmIx6w9o+RZLTuHxH35R/fAdyHmZmZ2bZbGREXjHUjHrD2oBxfdSnpr5ovR8RFdRmt00ln5n1N0k7AZ/NUgHPyJVsnASuBY0bYYRueRetUIuqoApFRVYjeqoSy78cK/L6U/CMHgIEKfG87qPJHbw1S7hrV6cOdnT9hqhVa8hvrAWuPKcRXHUU6sepWSQsjYkUto1XSDrVEAElzgTPyOh8DrgYOBA6MiOVdeRJmZmbWUAX+Fqr31lZspAJ/S9sobXN8VUSsjoj/ANZ0qlgzMzPrXRGxqhXb8RHW3uP4KjMzs15VvSOsLeEjrH0qIi6LiFcAnyDFV42YpPmSlklatn7d6vYUaGZmZi+haM2tajxg7T3bEl913Gg6iIgFETErImaNnzB1G0o0MzMzGzkPWHvP5vgqSRNI8VULiw0k7V/40fFVZmZmVVG9S7O2hOew9pi6+KpB4PKIWC7pAmBZRCwEzpL0NmA98DSFq1lJegjYAZgg6Tjg6IhY0ennYWZmZnUqOthsBQ9Ye9MmtrytNwJExHmFx+9my4lW43M7JB0FPAX8J7AO+HjTwarKn3Na9vrANbaKa2yBgfL/bzgwUO5sToBx2tjtEpoa7xrHrAKRwD3DA9YeM1wOa6HZNRHxxdx+LnAJMAd4EviTiHhU0oGko7R7dPQJmJmZWUOifwfJZf9b30ZvLDmsP4+IR/Py5cBkSRM7ULOZmZmNhOewWo9oVQ7r8cDPImJtg3XnA/MBJk6e1oKSzczMzIbmI6x9argcVkmvAT4L/Lch1nWslZmZWRc4h9V6xZhyWCXtCXwbODUiHmhLhWZmZrZt+nRKgAesvWebc1glTQOuB86NiB93qF4zMzOzYXkOa4/JOaxXAPeQTia8qUEO65ckHUn6G2stcEpe/TPAgcDXJAE8BhwREU8M2+dAuc9ZLHt9UIGoI6pRYyX+BC97bFT5f10YrMD3meMqEL01oPLXOFjFQ4Ht1qe7pAof7zYKOdbqNOBVpASA3SXNjIjz8mAV4LiImBQRk4FTgTPz8rOBiXn5fnn9VR19AmZmZtZYi+avVuBvvpfwgLX3jCXW6oWI2JCXT6Jv/44zMzOzMvGUgN4zplgrSUcAlwN7A6cUBrBmZmbWbX16KMlHWPvUULFWEXFLRLwGOAz4pKRJ9etKmi9pmaRl69et7lzRZmZmfc5TAqxXjCnWqiYi7gaeJ52EVf+Yc1jNzMysYzxg7T1jibXaV9K4fH9v0olbD3WiaDMzMxsB57BaL8hzTmuxVquBx2qxVpLm5mZfkrRG0ovA1cD5efnrgTskLQd+BSyOiCc7/BTMzMxsCP06JcAnXfWYulirR4Bba7FWhWbH1ZIC8iD2TOB7EXE1cLWk64AVwPJm/YUgBlv8JFqtArmSzjhtjSrsx7LXqLLnxAKqwP+2AxWosQoZp2WvUSWvr5d4wNp7NsdaAUiqxVqtqDUYKtYqtz8OeJB0dNbMzMzKoqJf57eCB6y9Z5tjrSRtR0oNOAo4p+2VmpmZ2ej06YC15F9OWbsMEWv1aeDvI+L54dYtxlptWOsDsWZmZtZePsLae7Yl1uoL+f4RwLslXQxMAzZJWhMRny+uEBELgAUA203fs0//1jMzM+ssUc0TplrBA9besznWijRQPRE4udhA0v4RcV/+cXOsVUS8odDm08Dz9YNVMzMz6yIPWK0XRMQGSWcBi4FB4PJarBWwLCIWAmdJehuwHngamNe9is3MzMyG19YBq6Q5wKWkgdOXI+KiusdPJ0UqbSRdVWl+RKyQtA9wNylLFGBpRJye13kd8FVgMrAI+EhEhKSdgG8C+5DC7k+IiKclKdfwduAF4LSI+Fne1jy2zN/8m4i4sof6WFTrQ9L7cx/fy23OB2bmPn4HPJq3807gQmATsAH4S5oSMVDu3KiyxwiBa2yVKtRY+pi1CsRaDVagxnHa1O0Smhqvjd0uoakqxIN1mqI/90nbPt5zHuhlwLGkwdFJkmbWNbsmIg6KiEOAi4FLCo89EBGH5NvpheVfAD4E7J9vc/Lyc4ElEbE/sCT/TO6/1nZ+Xp88MPwUad7m4cCnJE3v8z6WAH+QX48PAF/GzMzMyqFVV7mq4Ji3nccjNueBRsQ60sk97yw2GC4PtBFJLwN2iIilERHAVcBx+eF3Alfm+1fWLb8qkqXAtLydY4AbI2JVRDwN3AjM6ec+IuL53HZEr4eZmZl1Vr9e6aqdA9ZGeaB71DeSdKakB0hHWD9ceGhfST+X9G+SaicD7ZG302ibu0XEynz/MWC3JnUMt7xf+0DSuyT9EriedJTVzMzMrKu6PuNriDzQlcDvRcRrSeH210jaYRTbbPsB717tIyK+HRGvIh11vbDROsUc1vVrh41sNTMzs1bylICW25Y80NpX02sj4ql8/zbgAeCAvP6eQ2zz8fw1eG3qwBNN6hhueb/2sVlE/AjYT9KMBo8tiIhZETFr/MTt6h82MzOzNvGUgNbbnAcqaQIpD3RhsYGk/Qs/bs4DlbRLPmkLSfuRThj6Vf4a+zlJs/NZ86cC383rL2RLPNO8uuWnKpkNPJu3sxg4WtL0fJLS0cDifu5D0itzWyQdCkwEnsLMzMz6jqQ5ku6RdL+kc4dpd7ykkDSrXbW0LdZqjHmgbwQukLSeFLF0ekSsyo+dwZY4qBvyDeAi4FpJHwQeBk7IyxeRoqDuJ8VBvT/Xt0rShaSBNcAFvdIHcLikS4EpwHeA1XV9PAacL+l84EVSAgHA8cCH8lFXkaZmTATWYGZmZt3XoaOjhbSno0jnwdwqaWFErKhrtz3wEeCWttYTfZrn1avyG+xeCm8w4KTiG0zSDrWEBklzgTMiYo6kccDPgFMi4g5JOwPPRMSQYX3b7bRXHPzWj7TxGY3dhsllD76E9VPLX+O6Hcpf4/rty/95tm5aufM5NX1dt0toaped/7PbJTR1wPSXzLQqnSN2fLDbJTT1lin3NG/URSf88e9Yfue6jn04Tt15rzjwHR9tybZ+evXZt0XEkEdEJf0h8OmIOCb//EmAiPhMXbt/ICUUfRw4JyKWtaTAOl0/6cpabixxYkcDd0bEHbndU8MNVs3MzKxnNU17ytMH94qI69tdjC/N2nsavcGOqG8k6UxSAsME4C158QFASFoM7AJ8IyIubm+5ZmZmNmKt+yJphqTi0dAFEbFgpCtLGiBd8Om0llU0DA9Y+1REXAZcJulkUpzYPNL74fXAYaR5sksk3RYRS4rrSppPutoWE6ZM62jdZmZm/Uq09Az/J4ebEkDztKftgQOBm/P52rsDCyXNbce0AE8J6D3bHCdGOhr7o4h4MiJeIJ3odWj9Co61MjMz63nDpj1FxLMRMSMi9omIfYClQFsGq+ABay/a5jgxUqLDQZKm5BOw3gRsdTagmZmZdVFEa25Nu4kNQC3t6W7g2lraUz5hu6M8JaDHjCVOLCKelnQJadAbwKJOTKQ2MzOzkelk6H9ELCJ921pcdt4Qbd/czlo8YO1Nm9hy8bWN8JI32N1sOdFqfG6HpPcC55AGsgDnSLomIm4ftreSpx3FQMkLBKIK33WUfzdWYz+WvEaVvD4AVeAyPeNV7vgygPEqfwjMYMlf6wp8LPaMCnw02WgUgn6PBWYCJ0maWdfsmog4KCIOAS4mneVHRHwtIg7Jy08BHmw6WDUzM7POiBbeKsZHWHvP5hxWAEm1HNbNc1GHyWEtOol0QpaZmZmVRAUO3reFB6y9Zyw5rEV/Rt0FB8zMzMy6wVMC+lREXBYRrwA+Qcph3UzSEcALEXFXo3UlzZe0TNKy9Wuf70C1ZmZmBvTtlAAPWHvPWHJYa04Evj7UCs5hNTMz6w5Fa25V4wFr7xlLDmvtUmsn4PmrZmZmVhKew9pjcg7rFcA9pMSNmxrksH5J0pGkLwXWkhIBkDQe+B4wHbhe0lUR8ZmmfZY8NqoKUUdl34dQlf3Y7Qqai4FyH9qoQmRU2aOOAAYqUOMg5T97Z7Dk312r0/UFIwr970UesPaYHGt1GvAq0glXt0qaWZfDelwtKSBfreJM0kD1PcAzETFF0hRghaSvR8RDnXwOZmZm1lgF/hZqiwocj7BR2hxrFRHrSF/tb3W2/zCxVgFMzZdlnQysA4ptzczMzDrOR1h7z1hira4jDW5XAlOAj0bEqrZWa2ZmZiPnI6zWT4aItTqcdCnXlwP7AmdL2q9+3a1jrVZ3rGYzM7N+JpwSYL1jLLFWJwPfj4j1EfEE8GNgVv0KW8daTW1R2WZmZmaNecDae8YSa/Vr8vQASVOB2cAv216xmZmZNRfRulvFeA5rj8mxVmcBi4FB4PIGsVZnSXobsB54GpiXV78MuELSctI3D1dExJ2dfxZmZmbWSBW/zm8FD1h70ya2XHxtI0BdrNXdbDnRajxbpnCvA57P628iHa0dnsqffVn2+qAiNZY/KrYa3xmVfD9WIYe1Chmn4wc2druEpgZV/hzW8WX/fel2AX2kCh/vNgo5h/Uy4FhgJnCSpJl1za6JiIMi4hDgYuCSvPxDABFxEHAU8L/zla/MzMysDKJFt4rxYKT3jCWHdSZwU27zBPAMDU66MjMzs+5wSoD1ikY5rHvUN5J0pqQHSEdYP5wX3wHMlTRO0r7A69g6caC27pZYqzWOtTIzM7P28oC1Tw2Rw3o5aYC7DPgH4CfkObB1626JtZrkWCszM7OOCGBTtOZWMT7pqvdsSw7rFyAlDAAfrT0g6SfAvW2o0czMzLZF9caaLeEjrL1nm3NYJU3J+atIOgrYEBErOlO2mZmZNdOvc1h9hLXHNMlh3QScBOwmaS2wkq1zWHcFFudB6x7AnzTtj/JHMlUhjqns+xAqUmO3CxiJkr8fNVD+vTg4UP44poEKvBvHq/zRWxX42LEO8YC1B0XEImBR3eLzSV/vH0Wap3orcHLxCGpEPCRpFnB9bvNYZyo2MzOzEangVapawX+89I+mcVfZhcBngTWdLM7MzMya69cpAR6w9o+mcVeSDgX2iojrO1mYmZmZ2XA8JcAAyFe0ugQ4bQRt5wPzASZMmd7ewszMzCyp6FWqWsED1v7RLO5qe+BA4GZJALsDCyXNjYhlxQ1FxAJgAcDUnffq018dMzOzzhIgz2G1Hjds3FVEPBsRMyJin4jYB1gKvGSwamZmZtZpbR2wSpoj6R5J90s6d5h2x0uKfIY6knaW9ENJz0v6fF3b10n6Rd7mPyofDpS0k6QbJd2X/52elyu3u1/SnXmeZm1b83L7+yTN6+U+SF/31+Ku7gGmAt+R9CtJJ+X278zbvp10WdY/GNkrbWZmZh2xqUW3imnblABJg8BlFGKUJC2sD6KXtD3wEeCWwuI1wP8kfUV9YN2mvwB8KLdfBMwBbgDOBZZExEV5cHwu6bKjxwL759sRef0jJO0EfAqYRZoRcluu7+ke7uOGiDhA0sXAqkIfhwBfB5YACyMiJB0MXAt85aWv7tZKn89Z9vqg9NmcQDVqrMJrXfKc04GS1wcwUIFTnAdU/hHBYAVGLYPdLqCJbnwsekpA621zjFJErI6I/6AuWknSy4AdImJpRARwFXBcfvidwJX5/pV1y6+KZCkwLW/nGODGiFiVB3c3AnP6uY+IeD63hXQEtj9/K8zMzKxU2jlgbUeM0h55O422uVtErMz3HwN2a1LHcMv7tQ8kvUvSL0kXD/gAZmZmVg7RwlvFdO0LtEKM0tmt3nY+StjWl6NX+4iIb0fEq0hHXS9stI6k+ZKWSVq2Ye3qdpZnZmZmm0W60lUrbhXTzgHraGKUHgJmk2KUZjXZ5p5DbPPx/DV4berAE03qGG55v/axWUT8CNhP0owGjy2IiFkRMWvcxKn1D5uZmZm1VDsHrC2PUcpfYz8naXY+q/5U4Lv54YVA7Qz5eXXLT81n2c8Gns3bWQwcLWl6PhP/aGBxP/ch6ZWFtIJDgYnAU0O9HmZmZtZZ/Xpp1ralBETEBkm1GKVB4PKIWC7pAmBZRCwcbv181HUHYIKk44Cjc8LAGcBXgcmks+pvyKtcBFwr6YPAw8AJefki4O3A/cALwPtzfaskXUgaWANcEBGr8v1+7eN40qB4PfAi8GeFk7DMzMys2/r0v2V5PNJ7JM0BLiX9ofDliLhoiHbHA9cBh0XEMklHkQazE4B1wMcj4qbh+po6Y6+Y+ccfbWn9rbZhcrcraG791PJnRq3frtsVNLd+h/LH9GzccWO3SxjW5OkvdruEpvaa/ky3S2jqwGmPdruEpo7Y/oFul9DUmyeXez8e8/YnueOOdR37AN9h+z3i8Nee0ZJtLfn3v74tIoabhlkqvjRrjxlj/u2TwJ9ExKOSDiQdHd8q2cHMzMy6JKACEb9tUYWYbRudseTf/jwian/OLgcmS5rY7oLNzMxshJwSYD2iVfm3xwM/i4i1rS/RzMzMbOQ8JaDPFPJvTxumzWtIR1+PHuLx+cB8gAlTp7e+SDMzM2usegdHW8ID1t4zmvxbgN1J+bdz84lXewLfBk6NiIYz8iNiAbAA0klXrX8KZmZm1ogq+HV+K3hKQO/Z5vxbSdNIl2Q9NyJ+3I3izczMzOp5wNpjImIDcAVwD7AaeKyWfytpbrFtjrV6E/DqvOgzpKOvX5P0oqQHJe3awfLNzMxsOH160pWnBPSYHGt1GvAqtsRazYyI8+raFWOt7s6Lzwb+Il/04WXAHcAqmoiSR4hGBf4sq0KNlfjztuTvRQAGyv0fxUDJ6wNQBS7TM1CBGieo3JnAUP6PnY5/5ATgWCvrEWOJtXohH6EFmETfTu02MzMzSXMk3SPpfknnNnj8Y5JWSLpT0hJJe7erFg9Ye8+YYq0kHSFpOfAL4PTCANbMzMy6SASK1tya9rXlQkTHAjOBkyTNrGv2c2BWRBxMunLmxS1+ypt5wNpnCrFWZzd6PCJuiYjXAIcBn5Q0qcE25ktaJmnZhjWr21uwmZmZbdG5OaxNv7GNiB9GxAv5x6WkZKK28IC194wm1uohYDYp1mqr6wlHxN3A87ktdY8tiIhZETFr3KSpLS7fzMzMOmBG7eBTvs2ve7zpN7Z1Pgjc0Ooia3zSVe/ZHGtFGqieCJxcezAingVm1H6WdDNwTo612hf4TT7pam/SiVsPdbB2MzMzG07rzvB/MiJmNW/WnKT3AbNIyUNt4QFrj8mDzbOAxcAgcHkt1gpYFhELh1n99cC5ktaTzkM8IyKebH/VZmZm1lRnUwKafWMLgKS3AX8FvKmdl3P3gLU3bSK9rQPYCFAfawVb5bCS21wNXC3p94AVwCtH0lnZI5lioAJZRyXfhwBRgbijsr8XgdJHb1UhMmrcQPlzfcZXIjKqCvux3L/UJf91Hqthv7EFkPRa4EvAnIh4op3FeMDaYwpn9R3FlhzWhRGxoq5dMYe13iW0cR6KmZmZbZtOXZp1hN/Y/h2wHfCtfLn3X0fE3CE3OgYesPaezWf1AUiqndW3oq5dLYf148WFko4DHiRdJcvMzMzKpINXqYqIRcCiumXnFe6/rVO1lPtYu22Lbc5hlbQd8Ang/HYXaWZmZqPVokgrX5rVyq6Qw3pag4c/Dfx9RDyfD+0PtY35wHyACVOnt75IMzMzswIPWHvPaHJYAXYn5bDOBY4A3i3pYmAasEnSmoj4fLGDiFgALACYOmOv6v2ZZmZmVkVBJY+OtoIHrL1nm3NYgTcUln8aeL5+sGpmZmZdVP5wh7bwgLXH5LP6rgDuISVu3DRUDmsh1urVwDJJ+wB353V3J5189blh+1MFooQqkDsSrrElSv9eBFTyeLCBCsRaucbWGFT5Rz6DVfgAt47wgLXH5Fir00hXqarFWs2sz2Gti7W6u/DQAxFxSIfKNTMzs1HoVKxV2VTgeISN0uZYq4hYB9RirerVYq3WdLI4MzMzG4M+TQnwgLX3bHOsVbavpJ9L+jdJb2jwuJmZmVlHeUpAn2kSa7US+L2IeErS64DvSHpNRDxXt43NsVbjt3OslZmZWUcEsKl6R0dbwUdYe89oYq0eAmaTYq1mRcTaiHgKICJuAx4ADqjvICIWRMSsiJg1btLUNj0NMzMz21r/XjjAA9besznWStIEUqzV5mSAiHg2ImZExD4RsQ+wFJgbEcsk7ZJP2kLSfsD+wK86/xTMzMzMtvCUgB6TY63OAhYDg8DlQ8VaNfBG4AJJ60lJb6dHxKr2V21mZmYjUsGjo63gAWtv2kS+HgawEaA+1gq2ymElt/m/ku4DvgTsAPytpBsjYvgkgZIfp69CNmcVaiz76wxAyTNOgdLXWIX80HHa2O0SmqpCxul4yr8fB0r+waNu5MR6wGq9IH+lfxlwFFtyWBdGxIq6dsUc1tqyccD/AU6JiDsk7Qys71jxZmZmZg2U+08X2xZjyWE9GrgzIu4AiIinIqL8f4KbmZn1g1pKQCtuFeMBa+8ZSw7rAUBIWizpZ5L+e6MOJM2XtEzSsg1rVreydjMzMxtSQGxqza1iPCWgzzTJYR0HvB44DHgBWCLptohYUmwUEQuABQBTdtmren+mmZmZWaV4wNp7RpPDCrA7KYd1Lulo7I8i4kkASYuAQ4GtBqxmZmbWJX160pWnBPSebc5hJUVhHSRpSj4B603Aipd2YWZmZh3Xx3NY23qEVdIc4FJSHuiXI+KiIdodD1wHHJYD7N8LfLzQ5GDg0Ii4PV8y9KvAZGAR8JGICEk7Ad8E9gEeAk6IiKeVDiNeCryd9DVqimlZAAARcElEQVT3aRHxs9zvPOCvcx9/ExFX5uWV7gOo5bBOBsYD35X0FPC/ImJhXR87k466QsphnQE8lX/+5wbzXLem8kcylb0+gG4ko4yW92NrqOSxUQMD5Z/bVoXorfGO3mqJQVXgl9o6om3/BRXilY4FZgInSZrZoN1L4pUi4msRcUhEHAKcAjwYEbfnh78AfIh0Fab9gTl5+bnAkojYn/QV9rl5+bGFtvPz+uRB26eAI0hn1n9K0vRe6AP4fxFxALASOC6vs4otEVXFPi4u9L0E2DsiJuftvQ4zMzMrD1+ateXGEq9UdFJeF0kvA3aIiKUREcBVpAEZedtX5vtX1i2/KpKlwLS8nWOAGyNiVUQ8DdwIzOnnPiLi+dwWYCrpywczMzMrCw9YW24s8UpFfwZ8vbDNR4bY5m4RsTLffwzYrUkdwy3v1z6Q9C5JvwSuBz6AmZmZWZd1bVZaIV7p7GHaHAG8EBF3jWbb+ShhW/986NU+IuLbEfEq0lHXCxuts1UO64vOYTUzM+uMFh1d9RHWrYwmXukhYDYpXmlWoc2JbDm6WtvmnkNs8/H8NXht6sATTeoYbnm/9rFZRPwI2E/SjAaPLYiIWRExa9zkqfUPm5mZWTsEsGlTa24V084B61jilWpHYE8gz1/N66wEnpM0O581fyrw3fzwQmBevj+vbvmpSmYDz+btLAaOljQ9nwh1NLC4n/uQ9MrctjZdYyJbEgPMzMzMuqJtsVYRsUFSLV5pELg8IpZLugBYFhELh98CbwR+ExG/qlt+BlvioG7IN4CLgGslfRB4mDTYhRQZ9XbgflIc1PtzfaskXUgaWANcEBGreqEP4HBJlwJTgO8Aq4fo48PAy4C35uXvAT4haSLp77irCidhmZmZWbf16X/L8nikt+Q4sXuBo0gnWt0KnBQRK+rabU86sWoCcFbOvz2ZdJT7RElTSBcNeHNEPDRUf1N23SsOeM9H2/NkWmTD5PLn+G2Y0u0KmtuwXfk/K9bvUIHsyx3XN2/URdN3LP+89FdMf7LbJTR18Pa/bd6oy46cel+3S2jqv0wq9+/LkXN+y213rO3YfzI7jt81jtzp+JZs6/tPfPG2iJjVvGU5VCEK3EZnLHFiAUzNV7maDKwDnmtzvWZmZmbD8oC194wlTuw60vSBlcCvgc8VpjCYmZlZV7Xosqy+NKuVXSFO7LQGDx8ObAReDkwH/l3Sv9bPI5Y0n3S1LcZvN/0lGzEzM7M2CIio3hn+reAjrL1nLHFiJwPfj4j1EfEE8GPgJfNbHGtlZmZmneQBa+8ZS5zYr4G3AEiaShrM/rLTT8DMzMyG0KdTAjxg7TERsQGoxYndDVxbixOTNLfJ6pcB20laThr4XhERd7a3YjMzMxuxPr3Sleew9qZNpDP+gzQnlYg4r76RpOOBNxUWrQOez+tvYku267BioNyxUVGFP8vKvQsBiArUWIX9iMr9H8XgQPnnxw2094rVLTF+YEO3S2hqQFV4ravwS22d4AFrj8k5rJdRyGGVtHCIHNaPALcUFn8IICIOkrQrcIOkw6JfZ3ibmZmVSUQlL6vaClU49mSjM5Yc1pnATQD5pKtnaHDSlZmZmXVJn04J8IC194wlh/UOYK6kcZL2BV7H1okDZmZm1kWxaVNLblXjKQF9pkkO6+XAq4FlwMPAT8hzYOu24RxWMzMz6xgPWHvPaHJYAXYn5bDWoq0+Wmso6SfAvfUdRMQCYAHAlF33qt73CmZmZpVUza/zW8ED1t6zOYeVNFA9kXRBACDlsAIzaj9Luhk4JyKWSZoCKCJWSzoK2FB/spaZmZl1SVDJDNVW8BzWHtMkh/XTku6RdL+kcxusfgbwrKQ1wD8D/9K5ys3MzMwa8xHWHhQRi4BFdYvPJ329Xx939eZCmyeBL0bEWSPuTJQ/+7Ls9VGNrNgq1FiF13qg5DmsKnl9UP59CDBYgazYCS89RaF0BjXY7RKGpW586PRp0qQHrP1jc9wVgKRa3JW/8jczM6uAAMJTAqzHNY27yo6XdKek6yQ50srMzMy6zgNWK/oesE9EHAzcCFzZqJGk+ZKWSVq24cXVHS3QzMysb0WkKQGtuI2ApDnDnfsiaaKkb+bHb5G0T4uf8WYesPaPZnFXRMRTEbE2//hl0oUDXiIiFkTErIiYNW7y1LYUa2ZmZi8Vm6Ilt2YKl3o/lnQlzJMkzaxr9kHg6Yh4JfD3pCtotoUHrP1jc9yVpAmkuKuFxQaSXlb4cS4pZcDMzMz6z0gu9f5Otnwbex3wVuWQ91bzSVd9IiI2SKrFXQ0Cl9firoBlEbEQ+LCkucAGYBWNr4ZlZmZm3dK5lIBG574cMVSbPM54FtiZlDrUUoo+vWKCtYak35Eu49oqM2jDG73FXGNruMbWcI1jV/b6wDW2Sqtr3Dsidmnh9oYl6fsULv4zRpOANYWfF+QrWdb6ejcwJyL+a/75FOCIYvSlpLtym0fyzw/kNi1/H/gIq41Jq39RJS2LiFmt3GarucbWcI2t4RrHruz1gWtslSrUOJyImNPB7pqe+1Jo84ikccCOwFPtKMZzWM3MzMysXtNzX/LP8/L9dwM3RZu+uvcRVjMzMzPbygjPffkKcLWk+0nnvpzYrno8YLWyWdC8Sde5xtZwja3hGseu7PWBa2yVKtRYGo0u9R4R5xXurwHe04lafNKVmZmZmZWa57CamZmZWal5wGpdUabLvY2hxjdK+pmkDTn+o9P9ny7pQUlrJL0o6ZLCY5/M6z0h6WFJd0paIunURtvMk+pvycu/mSfYD/s6FPq4R9IxzeoueR+PSlon6be17TXqQ9LxkkLSXfV9FLZXazPqM5GHeV7F5V+TtELScknX5MdPk/ScpLX59pXCuq+T9Iu87j9KKdRb0k6SbpR0X/53el6u3O7+/L45tLCteXkfrZP0eGG/FPv4uaTb8+1+SesLffx5vn+fpB8U+nhvYf0fFNqcUahxaX7OtedRq3F5fqydfTws6QFJmyTNyvuhtv68Vu7rwmv9uKQnG/TxDknP59f5SUlvz8vfk2vfpG14742Gmn827a30eXOnpJsl7dnOehr0f7nSZ99dQzw+5HvcSiwifPOtozfS5O0HgP2ACcAdwMy6NmcAX8z3TwS+WcIa9wEOBq4C3t2F/qcV2vwp8Bzp8nkzc/uJwEnAr/L2zgCeb7RN4FrgxHz/i8CfD/c61PWxb65jcLi6S97HA8Ab8793AAc26GMW8CNSZMv/qO8j/7x9brMUmNWK17xu+UzgRWB2XmfX/O+ZwLPATsD0/JpPz4/9FJgNCLgBODYvvxg4N98/F/hsvv/23E55vVvy8p3ydh8EDsn378o1DdXHvwK35vvnA8/k7ZwAvJBrnU16X87Oj71AmhM3Pbc/P6//G+Dq3Mcy0hnMyst+3eY+Tsyvx83AW/Jzb8u+zq/1awv7d3ZdH0/k7Yn0Xns8L3818Pu5xlG999rw2fQtYF6+/xbg6g5/fr8ROBS4a4jHG77HfSv3zUdYrRtKdbm3ba0xIh6KiDuBdlx2ZCT76NW1NqRB3crc5p3ANyJibUR8Hbgvb281sL5+m3m/voW0nyHt9+Py/aFeh2IfDwL35z4a1l3yPn6St/ej3MdPSAPA+j7+gXSd7B2AJQ36ALgwtymGcY/UUK/55uWkq88tAv4IICKeyOseCPwmIlZFxNPAjcAcpcst7xARSyMiSH9cNdon9fvqqkiWAtPydo4hDU7ujYjbcx+/AN43TB//Bfhcvv9YrnkV8FbSfp5DuvDI+PzvMXn52/LzAHgs978GODz38QzwbL5/OLAut2lLHxHxDVLI+gTSAOfGNu3r3YFfA6/K2/1afq2LfYwHnsp9/CsQ+TnfHRH30H4j+WyaCdyU7/+wweNtlX+XVw3TZKj3uJWYB6zWDY0u97bHUG0iYgPp6NHOHamurv+sUY3d7n8PYDulK4tcDPxTXjbUuieSjozUL98ZeCbv5/q+hnodhupjqOVl7mNdYXuPAOtJR8+LfQySBhO3kC5dvLG+j/y14l4RcT3bZiTP9wDSHycfzl9f10LEpwN75683ryMdba+t+0iDbQLsFhEr8/3HgN1GUMcatt5XAbyiUR+S9s611v6AmJp/rvXxYGG7qwv3a8vJ7afmnx8u1DgFmFx7HvmxPdrYR+15TQR2HWb/jHVfP5tvteW1bRT7+AXwPkmPAB8n/UHaSSP5bLqD9K0PwLuA7SV18vO7mW5/vts28IDVrNrujohXAJ9g+KMYbyANLO7sSFU9RNIA8F7gx03aXAKc3eZyxpEGzt8mTff4J0nTgBXA5yLiYNLRuONHs9F8tK7VkTEnko7ob2zacuQ6EWvT1j5asK93B74aEXuSBqyH5PdfmZwDvEnSz4E3ka6G1Mr3gfWhsr3JrT+M5nJvqM2XexvCSGrsdv/FNt8ADsvL6tc9mDSY/Uu2PopQ2+ZTpK/ExtUt36qPutdhqPqGWl7mPiYUtrcn6SvXhwrLtgd+D5hLmtc4CViYT2yp9bE96Wv5myU9RPrauNZmpEbyfB/Jtf0mT5O4F9ifdJTt5bnNl0mvc23d4gkvxX3yeO1r0PxvbXrBcHVMYut9JdJR+0Z9nMiWr9ohHeFcW+hj38J2pxbu15aT26/OP+9dqPEF0lxegMfzY79tYx+157U2Pz7U/hnrvt4x32rLa9so9rE3aa42pKOxQeuuLT8STT+bIuLRiPjTiHgt8Fd52TOdK7Gpbn++2zbwgNW6oVSXextDjd3u/+laG9KAdFNusxA4UenM+HeQTpA5BvhBo23m/fpD0n6GtN+/m+8P9ToU+9iXNGj66VB1l7yPI/P23pD7OBL4QmHfvkia2/r6iNgH+B3w+YhYVusjIp6NiBkRsU9usxSYm9uM1FCv+eblwL+Q5ukulDSDNEXgV8DtwNH57POTSIOYxflr6Ockzc5zdk8dYp/U76tT85nUs0nzOFeSrnZzCPD7kv4AOBo4iDTPsr6P20jTFK4t9LE7QK7xpryfF5MGYOvzv4vz8iW5HcDuuf9JwE9zH9NIf5wo75+JuU1b+qjtB9L0kaW1fZ3XP7qF+7o2+L47b/e9pJOoin28CMzPffx5fl6/o3OafjZJmlE46vtJ4PIO1jcSQ73HrcyiBGd++dZ/N9JZmveSjs78VV52Aek/eUj/cXyLNFD4KbBfCWs8jHTEazXpaN3yDvd/KWnO3RrSf2KXFtpcnddbTRrY3p5vS+u3mdfZL+/n+/N+n9jsdSAdOXkAuId8NvRQdVegj5WkwcijhX19DemoS30fS0lnb2/VR91rdzPbcKb2MK/5+YXl/0GaAvA48He5zWdyrWtJRwY/WdjmrFzvA8Dn2XLBmJ1JJ4/dRzp5Z6e8XMBluf0vis8D+EDeR+tIRwlrNf4T6b1Y6+PTwEUN+viLvN/uz8trfZxSqHFJoc1ZhfVvIQ3kan3UalyRX9N29vFw3t9r8793FdZ/fyv3deE98ATpc+V+0gC39nv/HrYcSf4dcHRe/i7S51GtxsVd/Gx8d36u95KO+L/kd6TNn91fJ/1Or8/75IPA6cDpzd7jvpX35itdmZmZmVmpeUqAmZmZmZWaB6xmZmZmVmoesJqZmZlZqXnAamZmZmal5gGrmZmZmZWaB6xmZmZmVmoesJqZmZlZqXnAamZWMZIOk3SnpEmSpkpaLunAbtdlZtYuvnCAmVkFSfob0hW8JgOPRMRnulySmVnbeMBqZlZB+Trut5IuzXtkRGzscklmZm3jKQFmZtW0M7AdsD3pSKuZWc/yEVYzswqStBD4BrAv8LKIOKvLJZmZtc24bhdgZmajI+lUYH1EXCNpEPiJpLdExE3drs3MrB18hNXMzMzMSs1zWM3MzMys1DxgNTMzM7NS84DVzMzMzErNA1YzMzMzKzUPWM3MzMys1DxgNTMzM7NS84DVzMzMzErNA1YzMzMzK7X/D8b4R7ORpL5ZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "code", "metadata": { "id": "mcMsJ0QDDebo" }, "source": [ "" ], "execution_count": 5, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "I8TM6pzDDebo" }, "source": [ "" ], "execution_count": 5, "outputs": [] } ] }